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ABSTRACT 


This report presents the results of activities conducted over the period 1/2/85 - 
12/31/90, in which the study of forced convection boiling under reduced gravity was 
initiated. The study seeks to improve the understanding of the basic processes that 
constitute forced convection boiling by removing the buoyancy effects which may mask 
other phenomena. Specific objectives may also be expressed in terms of the following 
questions: 

(1) What effects, if any, will the removal of body forces to the lowest possible 
levels have on the forced convection boiling heat transfer processes in well-defined and 
meaningful circumstances? This includes those effects and processes associated with the 
nucleation or onset of boiling during the transient increase in heater surface temperature, as 
well as the heat transfer and vapor bubble behaviors with established or steady-state 
conditions. 

(2) If such effects are present, what are the boundaries of the relevant parameters 
such as heat flux, heater surface super-heat, fluid velocity, bulk subcooling, and 
geometric/orientation relationships within which such effects will be produced? 

A flow loop was designed and fabricated to permit operation at low velocities and 
various orientations in the test section. Flat heaters are used with flow parallel to the 
surface, and include both semi-transparent thin gold films on quartz and copper heaters. 
The gold films serve simultaneously as heaters and resistance thermometers, and permit 
visualization from behind the heater surface. Results are presented for both transient and 
quasi-steady nucleate boiling. The quasi-steady results do not include the full spectrum of 
orientations possible. 

The experimental data presented here are the results of the activities of the following 
Research Assistants: Dr. Jamie S. Ervin, Mr. Longhu Li, and Mr. Kevin M. Kirk. 
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1. INTRODUCTION 

This report presents the results of activities conducted over the period 1/2/85 - 
12/31/90, in which the study of forced convection boiling under reduced gravity was 
initiated. The study seeks to improve the understanding of the basic processes that 
constitute forced convection boiling by removing the buoyancy effects which may mask 
other phenomena. Specific objectives may also be expressed in terms of the following 
questions: 

(1) What effects, if any, will the removal of body forces to the lowest possible 
levels have on the forced convection boiling heat transfer processes in well-defined and 
meaningful circumstances? This includes those effects and processes associated with the 
nucleation or onset of boiling during the transient increase in heater surface temperature, as 
well as the heat transfer and vapor bubble behaviors with established or steady-state 
conditions. 

(2) If such effects are present, what are the boundaries of the relevant parameters 
such as heat flux, heater surface superheat, fluid velocity, bulk subcooling, and 
geometric/orientation relationships within which such effects will be produced? 

The potential effects, implied above, have their roots in observations of nucleate pool 
boiling under variable gravity perpendicular to the heating surface from high gravity to 
microgravity to negative gravity [16]. It had been observed that under high gravity 
conditions the nucleate boiling process is degraded; that is, for a give constant heat flux, the 
driving potential (heater surface superheat) is increased. For reduced and negative gravity 
conditions the nucleate boiling process is enhanced; that is, for a given constant heat flux, 
the driving potential (heater surface superheat) is decreased. These are illustrated 
schematically in Figure 1, where the observed temperature distributions with pool boiling 
in a saturated liquid are qualitatively presented. The variation of the buoyancy has an 
influence not only on the heater surface temperature, but on the boundary layer as well. 
The research, whose initial results are presented here, involves the determination of the 
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influence of an imposed velocity parallel to the heating surface on the bubble dynamics and 
on the resulting heater surface temperature and liquid temperature distribution. 

The enhancement of nucleate pool boiling with reduced gravity is believed to be due 
to the influence of buoyancy on the size and thickness of the microlayer trapped under a 
growing vapor bubble. Any residual influence of buoyancy on nucleate boiling in the 
presence of an imposed bulk liquid velocity, say parallel to the heater surface, depends on 
the extent to which the microlayer would be affected by the combination of buoyancy and 
forced convection. This is governed in turn by the various forces acting on the vapor 
bubble as the dynamic evaporation/condensation processes are taking place. Figure 2 
shows the various forces acting on a vapor bubble in nucleate boiling. With forced 
convection of the bulk liquid parallel to the heater surface two forces are acting in addition 
to those involved in pool boiling: the drag or liquid shear on and around the bubble as a 
result of the bulk liquid motion, and the lift generated by the liquid velocity change as it 
moves around the bubble. 

A total research program may be subdivided into three sequential phases, each 
intended to provide the base for the next phase: 

Phase A- This consists of testing in the laboratory with a flow loop at variable 
orientations between the flow direction and the gravity vector, using variables of flow rate, 
heat flux and subcooling with sizes, construction, and orientation of heating surfaces to 
serve as preliminary models for Phase B. It is expected that such testing will result in data 
that will provide guidance for the testing program at reduced gravity. 

Phase B: This would involve tests at reduced gravity using aircraft flying parabolic 
trajectories, with a portable flow loop and instrumentation resulting from the developments 
in Phase A. This would serve to determine the parameter boundaries for Phase C more 
accurately and economically. 

Phase C: This would involve orbital space flight testing in the shuttle or equivalent 
vehicle, with long term microgravity and well-defined parameters. The long terms 
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available will permit attaining a uniformity of experimental conditions between the various 
specific tests. It is foreseen that the test package used, resulting from Phases A and B, 
would be compact, self contained, and virtually completely automatic. 

The activities described below include only results from Phase A to date. Although a 
major em phas is is placed here on experimental measurements to observe the behavior of 
forced convection boiling under microgravity, appropriate analytical activities are an 
integral component of the total research program. 

1.1 General Background 

Phase changes with or without forced convection can provide high or low heat 
transfer rates, depending on the mode, i.e., nucleate or film boiling, evaporation, film or 
dropwise condensation. The mode is governed by a number of factors such as the degree 
of superheat or sub-cooling present in the liquid, vapor and container walls, the fluid/solid 
properties, body forces present, fluid velocities, and system geometry/orientation. 

Requirements for the proper functioning of equipment and personnel in the space 
environment of reduced gravity and vacuum introduce unique problems in temperature 
control, power generation, energy dissipation, the storage, transfer, control and 
condition in g of fluids, (including cryogenic liquids), and liquid-vapor separation. Boiling 
in microgravity is fundamentally different from boiling in earth gravity: the buoyancy force 
which induces liquid and vapor motion in boiling with earth gravity is effectively eliminated 
in microgravity. Temperature control in certain locations where internal heat generation 
tain* place, either as a result of dissipation or because of the nature of the process, may 
require that this energy be transported to other locations of the facility, where it can be 
eliminated by radiation to space. Fig. 3 illustrates two advantages in the use of phase 
change for the transport of energy in space; not only is the pumping power reduced by a 
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factor of 484 for the conditions assumed, but the mass of the fluid required is reduced by 
one-half. 

Certain effects which can be neglected at normal earth gravity, such as surface tension 
and vapor momentum, can become quite significant at microgravity conditions. Examples 
of applications in which these effects must be considered are: ullage control in storage 
containers; mechanisms acting in heat pipes; the effective transfer/flow of sanirntcd or near 
saturated liquids from one vessel to another. The latter is a particular problem with 
cryogenic liquids, where the transfer lines must be chilled, resulting in vapor production 
and two-phase flow. A phenomenon very similar to this in the mechanisms involved, and 
which gives rise to the fundamental study outlined in this proposal, is the process of forced 
convection boiling heat transfer under reduced gravity conditions. Applications in space 
stations are being considered [1], whether for temperature control or for vapor generation 
itself, as having distinct advantages over passive heat pipes in certain circumstances. 
Experiments of two-phase heat transfer conducted at earth gravity are generally either in the 
horizontal or vertical orientation. With the horizontal case separation of phases occurs due 
to gravity which, together with interfacial shear, can give rise to severe surging or 
chugging. This behavior may be quite different with significantly reduced gravity 
conditions. In the vertical orientation, either with upflow or downflow, the body forces 
accelerating or decelerating the vapor phase relative to the liquid will likewise produce 
behaviors quite different than in a microgravity environment 

Very small temperature differences, whether superheat or subcooling, which may 
normally be of little importance, can produce significant effects when the processes are 
diffusion lim it e d, as will be the case under microgravity conditions. Their influences must 
be understood, anticipated, and given appropriate consideration. Such small temperature 
differences can arise in large containers subjected to solar heating, for example, even with 
multi-layer insulation installed, where insulation penetrations are necessary for supports 
and fill lines. Additional heating can occur with connections to heated engine components. 
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Small temperature differences can also cancel the effectiveness of surface tension control 
devices such as screens. 

The effective application of forced convection boiling heat transfer in the microgravity 
environment of space, then, requires a well-grounded and cogent understanding of the 
mechanisms involved. Many correlations for a/g = 1 are presently available in the literature 
and are continuing to appear [e. g., 2-11]. This list is by no means exhaustive, and the 
mere existence of such a large number means that each has its inadequacies and limitations. 
Any attempt to extend these correlations to microgravity conditions, or to modify them 
using experimental results obtained without adequate consideration of mechanisms of a 
fundamental nature, as has been proposed [12], will provide results of limited utility. 

A discussion of the mechanisms in forced convection boiling heat transfer anticipated 
to be influenced by changes from earth gravity to microgravity will be presented below. A 
review of early works on pool boiling under reduced gravity is available [13], along with 
more recent data [14, 15, 16]. 

1 .2 Basic Mechanisms of Forced Convection Boiling 

a. Nucleation. The onset of boiling is inherently a transient process in the sense that 
once having begun, the dynamics of the boiling process will so change the situation that it 
can only be repeated by beginning anew. A number of studies of nucleation and the 
inception of boiling under non-forced convection circumstances with reduced gravity have 
been reported [17 - 20]. The condition at which nucleation takes place essentially depends 
on the microgeometry of the solid surface, the solid/fluid properties, the surface 
temperature of the solid, and the temperature distribution in the liquid. The latter two 
parameters in turn depend upon the imposed heat flux, whether saturated or subcooled 
conditions prevail, the velocity distribution in the vicinity of the heater surface, and the 
magnitude of the net buoyancy forces relative to momentum effects associated with the 
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velocity. The ratio of the latter forces constitute the Richardson Number, and is expected to 
be one of the parameters necessary to quantify the role that microgravity plays in forced 
convection boiling. The Richardson Number can also be expressed as the ratio of the 
Grashof Number to the square of the Reynolds Number (Gr/Re 2 ). For R,»l natural 

convection dominates, while for Rj«l natural convection effects can be neglected [113], 
and for Rjsl the full flow must be considered, and is referred to as mixed convection. In 

the absence of fluid turbulence this is amenable to computation. An example is included in 
Appendix A for the geometry used here. It is useful to define a corresponding "two-phase" 
Richardson Number by replacing the Grashof Number by the Archimedes Number, which 
is also a ratio of buoyancy to viscous effects, except buoyancy is now in terms of a finite 
density difference Ap: 

x f (i) 

A ’Two-Phase" Richardson Number thus provides a measure of the buoyancy versus flow 
forces in two-phase flow: 

RiW) - ^ (2) 

Results of research on nucleation at standard earth gravity have been reported for both 
pool boiling [21, 22] and forced convection [23, 24] conditions, while the thickness of the 
thermal layer at the initiation of nucleate pool boiling has been measured [25]. 

Once boiling has initiated and reached a steady condition, the nucleation site density 
becomes an important parameter in the description of the boiling. A reasonable amount of 
measurements of nucleation site density have been reported for pool boiling [26 - 29], but 
only one work is know for forced convection boiling [30]. Measurements of the tempera- 
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ture distribution in the boundary layer, which influences the nucleation site density, are 
available for pool boiling only [31, 32], and the interactions taking place between adjacent 
nucleating sites with pool boiling has been investigated [33]. The nucleation site density is 
also expected to be dependent on the departure size of the bubbles, and the factors which 
govern this will be considered below. 

b. Growth/Collapse. Once a particular nucleating site has become activated, the 
subsequent rate of growth and possible later departure and/or collapse are dependent on the 
transient temperature distribution in the vicinity of the bubble interface. The rate of growth 
affects the bubble frequency and together with the nucleating site density governs the 
relationship between heating surface superheat and the heat flux for pool boiling [32]. The 
rate of growth will be influenced by forced convection and reduced gravity only insofar as 
the temperature distribution is affected. Considerable work, both analytical and 
experimental, has been reported on the dynamics of vapor bubbles in the liquid bulk and 
near solid walls with pool boiling [e. g., 34 - 46]. The collapse of cavitation bubbles arc 
mechanistically the same as boiling bubbles [47 - 50], with large collapse rates associated 
with large subcoolings and large temperature gradients in the immediate vicinities of the 
bubbles. However, surface tension effects are neglected relative to the dynamic effects. 
With the slow velocities expected to be utilized with forced convection under the 
microgravity conditions of space, for energy conservation, it is not anticipated that the large 
liquid momentum associated with large collapse rates will be present, which result in 
cavitation damage. Instead, it is intuited that the growth and collapse rates will be relatively 
small, although still significant in their influence on the heat transfer rates, because the 
vapor formation arises from the relatively small liquid superheats. In this case the 
influences of surface tension may very well play a significant role in the heat transfer from 
the solid surface on which boiling is taking place. Again, a considerable literature exists on 
this effect with pool boiling [51 - 59], but none with the addition of forced convection. A 
factor in addition to surface tension which may become significant in the absence of body 
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forces is the momentum effect associated with the density changes of phase change [60, 
61]. This can influence the departure size of the vapor bubble as well as its subsequent 
trajectory, which will be discussed below. 

One further facet of vapor bubble nucleation and growth as influenced by surface 
tension should be mentioned here. The superheat that the liquid acquires in the boundary 
layer adjacent to the heater surface can be considerable, prior to nucleation. It is thus 
possible for the vapor formed initially to completely envelope the heater surface. With 
certain configurations such as small wires or cylinders it is possible that subsequent surface 
tension effects will maintain a stable "pseudo" film boiling process only because of the 
particular geometry used. It is expected that even if film boiling becomes suppressed to 
nucleate boiling on a small wire or cylinder, the thermophoretic effects and the resulting 
heat transfer will be quite different than with flat surfaces. Observations made that pool 
nucleate boiling is uninfluenced by changes from earth gravity to microgravity [62] are 
believed due to the large surface tension effects associated with the fine wire used, so that 
buoyancy is relatively unimportant in either case. The possibility of such effects, together 
with the fact that a flat surface provides a more well defined orientation for buoyancy and 
forced convection purposes provides the motivation for using a flat heating surface in the 
initial studies here. 

c. Departure Size and Trajectory. The size and trajectory of the vapor bubbles upon 
departure following growth in the vicinity of the walls will be important factors in 
establishing the flow pattern taking place in the bulk fluid stream, which can influence the 
subsequent heat transfer processes taking place as well as the pressure drop. For pool 
boiling the forces which play a role in the departure process are surface tension, buoyancy, 
inertia, pressure difference between the inside and outside of the bubble, and the 
thermophoretic forces resulting from surface tension gradients [63 - 68]. A possible source 
of error in assessing the departure of vapor bubbles from a solid surface has been pointed 
out recently [69], With forced convection, additional forces affecting the departure are 
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shear stresses and lift associated with circulation around the bubble, because of the velocity 
gradient in the flow field. 

The various forces acting are illustrated in Fig. 2. An analysis of the lift forces 
conducted for the case of potential flow is included here as Appendix B. 

Both analytical and experimental works have been reported [70 - 72], which include 
the possibility for sliding rather than departing, and with limiting effects taking place at 
very low velocities. Buoyancy was always present in these experimental works, of course, 
and a limitation exists in extending such measurements to behavior in microgravity 
conditions. 

The trajectory followed by a bubble following departure depends on the dynamics of 
the growth process, which will be influenced by the degree and distribution of subcooling 
in the flow stream, along with the fluid velocity gradient. The description of the motion is 
complicated by virtual mass effects [73 -76], and by interactions between bubbles and 
solid walls in the presence of temperature and velocity gradients [77-81]. It can be 
expected that the absence of buoyancy in microgravity will have a significant influence on 
these interactions. 

d. Pressure Drop. An extremely rich literature deals with the matter of pressure drop 
prediction in two-phase flow [e.g., 82 - 91]. In micfogravity conditions the pressure drop 
will be due solely to viscous and to acceleration effects associated with the quality changes 
with boiling. However, each of these will be dependent upon the size and spacial 
distributions of the vapor bubbles. As pointed out earlier, these depend upon the departure 
size and the subsequent growth and trajectories of the vapor bubbles which depend, in 
turn, on die temperature and velocity gradients. These are expected to be influenced by the 
removal of gravity forces, for the cases of low velocities to be used. 
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1 . 3 Effects of Gravity/Orientation in Forced Convection Boiling 


Little work has been done to investigate the effect of gravity in forced convection 
boiling, since the liquid momentum is normally assumed to be predominant over buoyancy 
effects. Some results available are described here. 

Bubble growth and motion of a single bubble with no gravity has been investigated in 
an isothermal and superheated fluid with a velocity field present [92], but the test time for 
these experiments were less than two seconds, and the study ignored the effect of a 
temperature distribution in the fluid as well as the interaction with other nucleation sites. 

The relative influence of velocity and buoyancy on the critical heat flux of liquid 
nitrogen was investigated [93]. The data presented was separated into buoyancy-dependent 
and buoyancy-independent zones. The inlet velocity required to prevent buoyancy from 
influencing the critical heat flux was found to be a function of the pressure and subcooling. 
This result is consistent with expectations, since the magnitude of the buoyancy is directly 
related to the volume of vapor present which is a function of the system pressure and 
subcooling. 

Transient slightly subcooled forced convection boiling experiments were conducted in 
a drop tower to simulate zero gravity [94]. The liquid velocities used were of the same 
order-of-magnitude as free convection velocities at earth gravity. It was observed that at 
zero gravity the vast majority of bubbles remained attached to the surface forming what was 
called a "bubble boundary layer." This phenomenon was peculiar to zero gravity, since 
bubbles always separated at earth gravity where the free convection velocities were nearly 
the same as the liquid velocities in the forced convection boiling. A correlation for the size 
of bubbles was obtained from a thermal equilibrium analysis and found to be a function of 
the saturation layer thickness. 

A few investigators have considered the effect of flow direction on forced convection 
boiling. The effect of surface orientation on bubble frequencies in nucleate pool boiling of 
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R-ll was investigated [95], and it was demonstrated that the bubble frequencies increase 
by increasing the orientation angle. Boiling nitrogen, was studied with upward and 
downflow [96], and the higher accumulation of void observed in downward flow 
suggested that different heat transfer coefficients may occur in upflow and downflow. 

A q uali tative comparison of upflow/downflow heat transfer of R-ll 3 was provided 
for a Reynold's number range from 1 x 10 4 to 5 x 10 4 [97]. These data show that the heat 
transfer coefficient for upflow is significantly greater than that for downflow in subcooled 
boiling. A corresponding but smaller difference also exists for saturated boiling. This 
difference between upflow and downflow is probably due to the fact that in upflow the 
buoyant and drag forces on a bubble prior to detaching from the surface are additive, while 
for downflow they are in opposite directions. Thus, it can be expected that bubbles detach 
from the surface at smaller diameters in upflow than for downflow, resulting in greater 
microconvection effect and enhanced heat transfer in upflow. In contrast to this work, Ref. 
[98] reported that for fully developed nucleate boiling with flow velocities of 0.2 and 0.8 
m/s the flow direction has a negligible effect on the heat transfer coefficient, even though 
the Reynold's numbers were lower than those of Ref. [97]. Photographs during boiling at 
low velocities show a higher frequency of vapor bubble formation in upflow compared to 
downflow, which was consistent with the observations of Ref. [95]. 

No obvious effect of flow direction on the heat transfer coefficient for fully developed 
nucleate boiling was observed in Ref. [99], where the entering liquid Reynold’s number 
was 1.4 x 10 4 . 
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2. HEAT TRANSFER MODELS AND CORRELATIONS 

Current theory on the mechanisms of heat transfer in nucleate boiling asserts that 
there is a component of heat transfer due to the rapid flow of heat through a liquid 
microlayer between the growing vapor bubble and the solid surface as well as a bubble- 
induced component due to the enhanced transient conduction which occurs as a result of the 
pumping action of the vapor bubble. In addition, with forced convection boiling there is 
also a contribution through single phase convection. There is considerable debate as to the 
importance of latent heat transport relative to the other mechanisms of heat transfer. 
Gunther and Kreith [100], using measurements from a photographic study of forced 
convection boiling, estimated that latent heat transfer accounted for only a small fraction (1- 
2%) of the total heat flow. Clark and Rohsenow [101] had similar findings. However, 
Bankoff and Mikesell [102] argue that if turbulent convective heat transport dominates in 
the heat flow of condensing bubble surfaces, latent heat may indeed be significant. 

Zuber [103] postulated that the mechanisms of heat transfer in nucleate boiling differ 
according to the different two-phase flow regimes. At low heat fluxes the vapor bubbles 
can be considered as isolated bubbles with no interference from either its predecessor or 
neighboring bubbles. In this regime, the heat transfer models based upon "bubble 
agitation" or "bubble pumping" give reasonable results. However, these models give 
incorrect results at moderate and high heat fluxes where bubble interference does occur. In 
this region of interference, the dominant heat transfer mechanism is by latent heat transport, 
according to Zuber. 

Experimental work by Bilicki [11] also supports the hypothesis that latent heat 
transport is significant in forced convection boiling. Bilicki argued that if heat transfer in 
nucleate boiling were enhanced by the bubble-induced turbulence, then the frictional 
pressure drop and the heat transfer coefficient should change simultaneously in accordance 
with the analogy between momentum and heat transfer. Yet, it was noted that the pressure 
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drop remained nearly constant during the onset of nucleate boiling, indicating that the 
bubbles do not act as a stirring device but as thermal sinks for the transport of latent heat. 

In a heat transfer model given by Del Valle and Kenning [104], an area of influence 
surrounding each nucleation site is considered, which is repeatedly quenched at the bubble 
frequency by liquid at the bulk temperature, and the heat transfer directly under the bubbles 
is modified by microlayer evaporation. In addition, the heater surface area between the 
areas of influence is assumed to be cooled by single phase convection. It was concluded 
that the enhanced transient conduction induced by bubble motion was by far the most 
important mechanism in heat transfer while microlayer evaporation accounted for only 2- 
3%, and single phase convection accounted for 5-10% of the total heat transfer. However, 
there is doubt as to whether the transient conduction model used accurately describes the 
quenching process, since the model assumes that the bubble-induced convection is 
sufficiently strong to produce instantaneous replacement of liquid in the quenching process. 

It has been assumed that the fully developed region of the boiling curve in forced 
convection boiling coincides with the extrapolation of the pool boiling curve. This 
assumption appears reasonable, at least for low velocities. However, experiments carried 
out by Bergles and Roshenow [23] demonstrate that the boiling curve for forced convection 
boiling is not a simple extrapolation of the pool boiling curve. Their forced convection 
boiling data does merge into an asymptote at large values of superheat, indicating the 
invariance of heat flux to velocity and subcooling in fully developed forced convection 
boiling. However, the slope and intercept of the asymptote differs from their pool boiling 
curve. Lemmert and Chawla [105] also noted that there was no significant effect of flow 
velocity or subcooling on the boiling heat transfer coefficient in fully developed forced 
convection boiling. 

Owing to the complexity of forced convection boiling, theoretical analyses cannot 
provide a general equation for boiling heat transfer coefficients for different substances and 
conditions. Hence, heat transfer calculations require use of empirical correlations, where 
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the correlations will have a relatively large uncertainty and can be used only for restricted 
cases. In the calculation of heat transfer coefficients, information on the operating 
conditions, the fluid properties, and the geometry are usually necessary. 

Chen [106] suggested that the heat transfer coefficient for flow boiling can be 
expressed as the sum of the heat transfer coefficients of forced convection, hf c , and pool 
boiling, hb. 


h = F(hf c ) + S(hb) (3) 

where the forced convection heat transfer is intensified by the factor F (F>1) an^ the 
boiling heat transfer is suppressed by the factor S (S<1). 

In addition to knowledge of the convective heat flux, qf C , and nucleate boiling heat 
flux, qb, some correlations also require knowledge of the conditions required for the 
inception of nucleation. Bjorge et al. [107] proposed: 


q = Vqfc 2 + (qb - Qbi) 2 (4) 

for subcooled boiling where qbi is the flux on the curve at a predicted superheat for the 
inception of nucleation, which does not depend on the available cavity sizes. A different 
superposition scheme proposed in this work obtains better agreement for subcooled 
boiling: 


q “ ' \j qfc 2 + qbi 2 ( 1 * { AT sat it/ATsai J 3)2 (5) 

where the recommended equations for the convective heat flux and boiling heat flux are 
given in Refs. [108, 109], respectively. The incipient boiling criterion for subcooled 
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conditions is derived by Bjorge [110]. The only empirically determined coefficient needed 
in the correlation is in the calculation of the nucleate boiling heat flux. 

Other correlations by Shah [1 1 1] and by Gungor and Winterton [1 12] consist of only 
the forced convection term, where the boiling effect is included in enhancement factors 
which are functions of the boiling number, Bo. 
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1 EXPERIMENT AL APPARATUS 


A small scale forced convection flow loop was designed and fabricated, occupying a 
total volume of 97 x 75 x 61 cm (38 x 30 x 24 inches), for use with R- 113 over a range of 
temperatures from 25°C (72°F) ambient to 60°C (140°F) with corresponding saturation 
pressure variations. An early version of the loop is shown in Fig. 4. As will be described 
in more detail below, the geometry of the conduit leading to the straightener section and of 
the preheater and condenser/cooler assemblies were changed. A schematic representation is 
given in Fig. 5, and provision is made to permit rotation of the assembly through almost 
360° relative to the gravity field vector while under continuous operation. To compensate 
for changes such as hydrostatic pressure taking place during rotation, the system pressure 
and temperature at the entrance to the test section are automatically controlled, as is the flow 
rate. To prevent cavitation at the pump inlet while operating near the saturated liquid state 
in the test section, heat exchangers are included for subcooling the liquid prior to pumping 
and flow measurement, followed by heating. The 12VDC centrifugal pump is capable of 
control over a 10:1 volume flow rate by the pump speed, using the output of a propeller- 
type flowmeter with a microprocessor to control the DC voltage. The outlet of the pump 
leads to a "T" section, one branch connecting to the pressure control and filling systems 
while the other branch leads to the heating system. 

The preheaters in the loop raise the temperature of the R-l 13 to the desired operating 
temperature. Referring to Fig. 5, subcooled R-l 13 leaving the condenser/cooler is pumped 
into preheater #1, a counter-flow heat exchanger. The heat exchanger raises the R-l 13 
temperature to within about 4°C of the set point temperature. Preheater #2 then heats the R- 
1 13 to the desired temperature level at the inlet of the test section, as indicated by a 
resistance thermometer at that location. 

Preheater #1 consists of a one pass multiple tube heat exchanger, with 35 stainless 
steel tubes, 0.953 cm OD x 0.699 cm ID x 33.02 cm length (3/8 in OD x 0.275 in ID x 13 
inches). Preheater #2 is identical to each of the three heat exchangers used as the 
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condenser/cooler, with these four (4) heat exchangers being similar to Preheater #1 except 
that each contains only 1 3 tubes instead of 35. 

A 2 kw a.c. heating unit in series with an Emerson (Model SA55CXJAR-4814) 1/3 
hp pump provides hot water to the shell side of the heat exchangers with R1 13 flowing 
within the tubes. A temperature controller maintains the desired hot water temperature by 
sensing the temperature of the water exiting the heater, and the R-l 13 outlet temperature 
from preheater #2 is controlled manually by a flow valve in the series with the pump. It is 
planned that this will also be automated in the future. 

The test section is shown in Fig. 6 and consists of a rectangular flow channel 4-1/8" 
wide, 14" long, with four different possible heights (1/8", 1/4", 1//2", 1"). Only test 
sections with heights of 1/8", 1/2" and 1" have been fabricated to date, which permit 
varying the bulk flow velocity by a total factor of over thirty (30) with the existing pump. 
However, the test section with a height of 1/8" results in difficulties in observation from the 
side, and its utility may be limited. The maximum Reynold's possible at present with R- 
1 13 in the 1" test section is about 5000. 

For this initial basic study, the rectangular flow channel in Fig. 6 provides for the use 
of flat heater surfaces, eliminating complications of surface tension effects associated with 
curved surfaces, and also provides a more well-defined flow field in the vicinity of the 
heater surface than would be possible were tubing or cylinders used. Additionally, the 
orientation between the surface and the gravity vector is more well-defined regardless of 
whether the gravity field is a residual one in space or earth gravity. 

The configuration shown in Fig. 6 permits the simultaneous use of three (3) pairs of 
identical heaters if the study of upstream or downstream interactions between boiling 
surfaces is desired, or the use of three (3) different pairs of heaters, either simultaneously 
or independently, without requiring the draining, disassembly and refilling of the loop 
system to change surfaces. This has been found to be a time consuming process. The use 
of heaters in opposing pairs within the channel permits simultaneous operation with 
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opposite orientation relative to the gravity vector and, where desired in the future, the study 
of the interaction between the boiling boundary layers. Further, the availability of three 
opposed pairs of openings in the test section permits the insertion of diagnostic devices into 
the flow stream if desired, such as local fluid velocity probes (Hot wire, LDV) or 
temperature traverses. The duct is a welded assembly of 1/2" al umin um plates with the 
various quartz windows and heater substrates held in place by bolted flanges with 
appropriate gaskets. A 10" long upstream section provides for smoothing and transition of 
the duct dimensions from the preceding 90° turn in the loop. 

Measurements of the velocity profiles in both the 1” wide and 1/8" wide test sections 
were made, using the hydrogen bubble technique in water. At the relatively low velocities 
used, in Reynold's number ranges from 300 to 1700, a maximum pump efficiency of 10% 
was measured, which is typical for these small pump sizes. Additionally, it was found that 
the velocity profile in the test section was quite unsymmetrical, owing to the attempts to 
make the loop compact in length. After a number of experimental trials involving relocation 
of the flowmeter, redesign of the flow turning section, which included the installation of 
flow turning vanes preceded by a large cross section flow duct filled with small chrome 
plated metal spheres to break up the liquid jet issuing from the tubing used, an acceptably 
symmetrical velocity distribution was obtained. 

The heated surface itself is rectangular, 19.1 mm x 38.1 mm (3/4" x 1-1/2"), in size. 
This size is large relative to the maximum anticipated bubble sizes, again to eliminate any 
geometrical dependency, and is small enough in an absolute sense so as to keep electrical 
heater power requirements within a manageable level. A representation of the flat heater 
surface relative to the flow is shown in Fig. 7. The heated surface is mounted on a circular 
substrate which can be rotated in its own plane, so that the heated length or aspect ratio in 
the flow direction can be changed conveniently by a factor or two. This will have the effect 
of changing the thermal boundary layer thickness for given levels of flow velocity and heat 
flux. Two different types of test surfaces are presently used, mounted as opposing pairs to 
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permit simultaneous or consecutive operation at a/g = +1 and a/g = -1. The first pair of 
heaters in the flow direction of Fig. 6 consists of semi-transparent gold coatings (400 A) 
on quartz to provide both steady and transient nucleate boiling, using the surface coating 
simultaneously as a resistance heater and as a resistance thermometer, and permit 
simultane ous high speed photography of the process from the side and through the heating 
surface as illustrated in Fig. 8. This permits obtaining data on the departure size and 
trajectory of the bubbles, along with the nucleation site density and frequency of bubble 
departures. A sketch of this heater is given in Fig. 9, and shows the means by which the 
current carrying and the potential lead electrical connections are carried through the quartz 
surface without introducing any impediments to the fluid flow. 

By passing a D.C. current through the thin film and measuring the voltage drop and 
current with sufficient accuracy the heat flux input and the mean instantaneous surface 
temperature can be measured, following appropriate calibration of the electrical resistance 
versus temperature. Accuracies of ±1°F have been attained with reasonable precautions, 
and its reliability and durability have been thoroughly tested. 

The heater can be used in several ways, referring to Fig. 10. If the fluid remains 
motionless, the temperature distribution in both the quartz substrate and the fluid can be 
computed from the measured power input, using classical techniques, and the measured 
surface temperature can be compared with the computed surface temperature. This is 
demonstrated in Figures 11 and 12 for a constant and uniform heat flux input in pool 
boiling, with the measurements following the 1-Dimensional computations up to the onset 
of natural convection at a/g = 1 in Fig. 1 1 and up to the onset of boiling at a/g = 10* 5 in 
Fig. 12. Once the fluid has been set in motion, whether by natural or forced convection or 
by boiling, or by both, the transient measurements of the thin gold film temperature and the 
power input as in the later time intervals in Figures 1 1 and 12 permit the computation of the 
mean heat flux to or from the substrate, and hence to the fluid. A limitation of the 
temperature measurement with the heater surface in Fig. 9 is that only the integrated mean 
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surface temperature is measured. With forced convection over the surface and a uniform 
imposed heat flux it can be anticipated that a temperature variation will arise in the flow 
direction as the thermal boundary layer develops. 

Under microgravity conditions and in the absence of boiling and farced convection in 
the fluid, the division of energy between the substrate and fluid can be computed using the 
well-known solution for transient conduction heat transfer in two semi-infinite solids 
having a plane heat source at their interface. For constant properties and constant plane 
heat generation rate at the interface, the division of heat between the two materials remains 
constant. 

The equations describing the transient interfacial temperature and division of heat 
transfer rate for a uniform constant heat generation rate at the interface are presented in Fig. 
13, along with the properties for quartz (Fused polycrystalline), pyrex, BK-7, R-l 13 and 
the resulting division of heat transfer rates between the R-l 13 and each of the substrates. 
For steady-state conditions the fraction of heat input transferred to the R-l 13 with forced 
convection and/or boiling can be determined once the steady heat loss through the substrate 
is known as a function of the interfacial and the surrounding substrate holder temperatures. 
This is obtained by calibration. 

Fig. 14 shows a copper heater surface with the same size as the gold film heater of 
Fig. 9, and is intended to provide a metallic substrate to the fluid more representative of 
engineering-type surfaces. Because of the large heat capacity this can be used only for 
quasi-steady operation. It is gold plated to provide the same heater-fluid surface energy 
relationship as with the gold-coated quartz surface. It also has the same external 
configuration as the quartz heater, and thus is interchangeable in the test section. A film 
heater is compressed against the underside of the copper body, which is then encapsulated 
within a stainless steel housing. A copper foil 0.001 inch thick is then soldered across the 
entire upper machined surface to eliminate the crevices which otherwise serve as false 
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nucleating sites. Heat losses through the underside are obtained by calibration, as a 
function of the heater surface temperature. 

A Heise Model 623 pressure transducer is used for the measurement and control of 
pressure at the inlet to the test section. This has a maximum pressure of 50 psi and a 
sensitivity to 0.0025 psi. An uncertainty in pressure of .025 psi corresponds to an 
uncertainty of approximately 0. 1°F in the saturation temperature for R-l 13. Temperature 
control at the inlet to the test section is maintained to within 0.1°F of the set point, which 
means that pressure control must be to within .025 psi to remain within 0.1 °F of the 
saturation temperature. 

For pressure calibration, a Ruska Pressure Calibration system was used. The air 
piston portion of this system is capable of resolving from 0.0001 psi at a level of 2 psi, to 
0.005 psi at a level of 600 psi. This was used to calibrate both a precision Heise Bourdon 
tube gage and the pressure transducer. The Bourdon tube pressure gage is used for in- 
place periodic check calibrations of the pressure transducer, using the filling connection to 
the test vessel. 

A 7" diameter stainless bellows is installed in an appropriate housing for the pressure 
control within the test section, and is shown in Fig. 15. Both visual and electrical 
indications of the bellows position are provided. This, together with the two pneumatic 
solenoid valves and a modulating proportional valve constitute the mechanical components 
of the pressure control system, which controls the steady state pressure to within +.025 
psi. The system is shown in Fig. 16. 

The flow control system illustrated in Fig. 17 serves to maintain the flow rate 
constant, as the system flow resistance varies due to changes in void fraction with boiling, 
or as the voltage output of the power supply changes. The turbine flowmeter, 
manufactured by the Halliburton Co., is rated for 0.3-3.0 GPM and has been calibrated 
over the full range, using water. 
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A schematic of the loop condenser/cooler control system, used to automatically 
maintain the desired degree of subcooling at the pump inlet, is given in Fig. 18. The 
condenser system also condenses any R*1 13 vapor which may exit the test section, and 
consists of two heat exchangers and a secondary pump. Heat exchanger #1 consists of 
three identical shell and tube heat exchangers in series, with each well baffled shell 5.7 cm 
OD (2.25 inches) containing 13 tubes 0.95 cm OD x 0.70 cm ID x 33.02 cm in length (3/8 
in OD x 0.275 in ID x 13 in length), and transfers heat from the R-l 13 in the main flow 
loop to water circulating through heat exchanger #1. Heat exchanger #2, which consists of 
a plexiglass shell and copper tubing formed into a spiral, transfers heat from heat exchanger 
#1 to cooling water supplied from the building. 

A degassing unit and storage tanks for the R-l 13, including a pressure gage on each, 
have been fabricated from 304 s.s. A typical assembly is shown in Fig. 19. The concept 
involved is a combination of distillation at room temperature, leaving behind high boiling 
point components such as oils and solids, and freezing of the R-l 13 on the fins using 
Liquid Nitrogen within the inner vessel. By maintaining a sufficiently low pressure, the air 
components, except for water vapor, remain in the gaseous state and are removed by the 
vacuum pump. A molecular sieve is installed prior to the freezing vessel to absorb the 
water vapor. 

A filling system, mounted on a portable can, has been fabricated for use with R-l 13, 
and consists of (1) a vacuum pump, including a trap, for removing the air prior to filling, 
(2) a flexible connection to the R-l 13 source tank, mounted on the can, (3) a Heise 
compound Bourdon tube-type pressure gage covering the range 15 psig to 25 psig, with 
minor divisions of 0.05 psi, easily readable to .025 psi, (4) a flexible connection to the test 
vessel, and (5) associated valves and fittings. 

Fig. 20 demonstrates the variety of parameters possible with the flow loop operating 
in earth gravity alone. Only a relatively few test results have been obtained to date, to be 
presented in the next section. Fig. 21 presents the definition of the orientation angle 
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between the flow direction parallel to the heater surface and the buoyancy vector, used in 
these results. For an angle -0 = 0°, the flow is downward, opposite to the net buoyancy 
acting upward. 0 = 90° means that buoyancy is perpendicular and away from the heating 
surface. 0 = 180° means buoyancy is in the same direction as the flow, while 0 = 270° 
means buoyancy is acting perpendicular into the heating surface keeping any heated liquid 
or vapor formed in the vicinity of the heating surface. 
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L experiment al procedi ires 

Because of the high solubility of air in R-l 13, precautions were taken in the filling 
process to insure that no contact was permitted between the degassed R-l 13 and air. A 
filling cart held a vacuum pump, a valve system and associated piping, and a pressure 
gage. The vacuum pump was used for the evacuation of the test loop and the connecting 
line. The piping connected the test loop to the tank holding the degassed R-l 13, and the 
valve system permitted the evacuation of the loop and piping without opening the valve of 
the R-l 13 container. These were evacuated with a vacuum pump for eight hours to remove 
air and water vapor. The pressure gage provided an approximate vacuum indication before 
the R-l 13 was allowed to flow from the storage container. In addition, this filling cart 
supported the elevated storage tank of the degassed R-l 13 which provided the hydrostatic 
pressure for filling. 

The storage container was irradiated with heat lamps to maintain a positive pressure 
relative to the atmosphere. When the pressure gage mounted on the storage container 
indicated a gage pressure of 34.47 kPa (5 psig), the valve between the storage container 
and the stainless steel lines was opened and R-l 13 flowed into the loop. The system was 
flushed with R-l 13 vapor several times before filling with the li quid , to further remove any 
possible residual gases. To avoid contamination of the degassed R-l 13 by atmospheric air, 
the pressure in the loop was always maintained above the ambient, either by the bellows 
actuated pressure system or by heating. 

For each test, a single point calibration of the thin film heater surface was performed 
before heating the loop to the desired operating temperature. As discussed earlier, the 
linear temperature-resistance relationship changed with time, with observable changes 
occurring after about a week. However, the slope of the temperature-resistance 
relationship remained constant, and hence, only a single calibration point was necessary for 
the determination of the new temperature-resistance relationship. Prior to the heating of the 
loop, the surface temperature of the heater was known with certainty, since the entire loop 
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was at the known ambient temperature, determined with calibrated thermocouples. This 
known surface temperature was determined simultaneously with the heater surface 
resistance for the calibration point data. 

The forced convection boiling loop was used to examine the effects of an externally 
forced velocity on incipient boiling with transient heating and various heater orientations, as 
well as on steady-state boiling. Consequently, both the gold film heaters with an imposed 
step increase in heat flux and the copper heaters with steady power input were used. The 
bulk liquid velocities selected were sufficiently low so as to not mask buoyant effects. 
Velocities of 1.7 and 4.5 cm/s produced corresponding Re Numbers of 508 and 1343 in 
the test section. For the non-boiling conditions, the maximum Richardson Number is 
estimated to be 60, computed from Equation (2) so that buoyancy effects may be expected 
to dominate. 

The experiments were initiated after the forced convection boiling loop had reached 
steady-state conditions of temperature, pressure, and flow rate. All tests were conducted at 
the nominal liquid pressure of 120.66 kPa (17.5 psia), with fluctuations of less than 0.7 
kPa (0.1 psia). The bulk liquid temperature was maintained either at 40°C (104°F), which 
resulted in a nominal subcooling of 12°C (22°F), or at 46°C (115°F), which resulted in a 
nominal subcooling of 7°C (13°F) for the imposed system pressure. Fifteen hours were 
required to reach steady-state conditions at the low subcooling level; four to five hours 
were necessary for the high subcooling level. 

Two orientations of the test section were used for the incipient boiling tests here: 
vertical and horizontal. The flow loop was rotated such that the heater surfaces were 
parallel to the earth gravity field (in the vertical orientation), and the buoyancy force was in 
the direction of fluid flow. The second test surface orientation was obtained by rotating the 
loop such that earth gravity was perpendicular to the heating surfaces. As a result, one of 
the surfaces was horizontal up, and the opposing surface was horizontal down, with 
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buoyancy aiding forced convection or hindering it, respectively. For tests conducted with 
steady-state boiling, a number of other orientations were used in addition to these. 


4.1 Spatial Mean Gold Film Surface Temperature 


The spatial mean thin gold film heater surface temperature, T w , is determined from 

the measured test surface resistance, R w , the single point calibration resistance, Rc, with 
its corresponding calibration temperature, Tc, and the slope of the calibration curve, in 

the following equation: 


Tw “ Tc + dR (V R c) < 6 ) 

J «T* 

For the constant heat flux tests, R w , Rc, Tc, and ^ were determined with 

representative values and uncertainties of 2.9985 ± 0.0008 ft, 2.6266 ± 0.0002 ft, 20.00 
± 0.01°C, and 214.83 ± 0.23°C/ft, respectively, as given in Appendix C.l. As a result of 
the uncertainties in these quantities, the heater surface temperatures were determined with 
an uncertainty of ± 1.0°C, as described in Appendix C-l. 

4.2 Metal Heater Spatial Mean Surface Temperature 

The metal heater surface temperature was determined from the chromel-constantan 
thermocouple imbedded in the copper heating block near the heater surface, as shown in 
Figure 14. The error in taking the temperature at the location of the thermocouple to be 
equal to the surface temperature was estimated to be less than 0.002°C, by ass umin g one 
dimensional conduction in the copper block and negligible contact resistance between the 
soldered copper foil and the copper block. An error of this magnitude is less than the 
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uncertainty of ± 0.05°C associated with the thermocouple temperature measurement and, as 
a result, was neglected. 

4.3 Bulk Liquid Temperature 

The liquid bulk temperature was measured by calibrated chromel-constantan 
thermocouples, with an uncertainty of ± 0.05°C (± 0.1°F). A discussion of the estimate of 
this uncertainty is in Appendix C.2. 

4.4 System Pressure 

The liquid pressure at the inlet to the test section was measured by calibrated 
transducers. The uncertainty of the pressure measurement was less than the desired 
uncertainty of ± 0.172 kPa (± 0.025 psi). The details of the estimate of the uncertainty are 
given in Appendix C.3. 
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1 EXPERIMENT AL RESULTS 


Experimental results obtained during the study period are presented below. The 
transient heating results are presented first, followed by the quasi-steady heating results. 

5.1 Transient Boiling 

The transient boiling experiments were performed using a step in heater surface heat 
flux. A step change in heat flux is the most elementary form of imposed heat flux variation 
possible, as all other time varying imposed surface heat fluxes can be constructed from a 
series of steps in heat flux. Since the resistance of the gold film heater surface is 
computed from the measured current and voltage drop across the entire heater surface, the 
resulting calculated surface temperature is a spatially averaged temperature 

The transient heater surface temperature for a representative test with pool boiling was 
shown in Fig. 1 1, along with a tabulation of the test conditions, and applies to the case 
where the heater surface is facing ^horizontal up in earth gravity. The onset of natural 
convection appears as an irregularity in the temperature-time plot. The time from the 
energization of the heater to the onset of natural convection is designated as tnc and is 
identified by the departure of the heater surface temperature from the one dimensional semi- 
infinite media transient conduction solution, and also by the observed onset of fluid motion 
as characterized by a wave-like disturbance recorded photographically. 

The next significant event in Fig. 1 1 following the onset of natural convection is 
incipient boiling, as indicated. Incipient boiling, or the onset of boiling, is defined as the 
appearance of the first vapor bubble on the heating surface. In some tests, this incipient 
boiling takes place at the point when the mean heater surface temperature reaches a 
maximum, while in other tests, the onset of boiling occurs prior to this point. For tests in 
which the latter took place, the level of the ensuing maximum heater surface temperature 
depended on the manner in which the boiling propagated across the heated surface. If 
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incipient boiling occurred as an almost explosive event over the entire heater surface, then 
the temperature associated with boiling incipience was the maximum surface temperature. 
If incipient boiling occurred at a location near a comer of the heater, for example, the heater 
would become cooled locally while the remainder of the heater surface continued to rise in 
temperature. The heater surface temperature measurement is a spatially averaged quantity, 
as described earlier, and the measurement would continue to rise until the boiling spread 
sufficiently to produce a subsequent decrease in the mean value. This boiling propagation 
will be defined in more precise terms later. 

Following the onset of boiling, a quasi-steady boiling region is noted in Fig. 11. In 
this domain of the temperature-time plot, boiling has spread across the entire heating 
surface, and the quasi- steady boiling temperature level is less than the maximum heater 
surface temperature, but above the saturation temperature for the liquid at the system 
pressure, as expected. 

Figures 22 - 24 present the measurements of the transient heater surface temperature, 
which is a spatial average, for experiments conducted with forced convection. The 
independent test variables included three heat flux levels, two subcoolings, and two flow 
velocities, in addition to the orientations described below. High speed video recordings of 
the incipient boiling across and through the heating surface were made. With both surfaces 
positioned opposite each other in the forced convection boiling loop, surface Q-13 was 
heated in the horizontal down position and surface Q-16 was heated in the horizontal up 
position. In the vertical flow experiments, where the forced convection boiling loop was 
rotated 90°, the R-113 flowed upward such that buoyancy assisted the externally forced 
flow. The data are tabulated in Appendix D. 

To illustrate the effect of orientation on the heater surface temperature, heater surface 
temperatures resulting from a nominal q" T of 4 W/cm^ and with similar initial conditions 

were used in Figs. 22 - 24. In Fig. 22, for the heater surface in the horizontal up 
orientation, the heater surface temperature decreased to a lower quasi-steady boiling level 
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once boiling had propagated over the entire surface, as also occurred in the case of pool 
boiling in Fig. 1 1 with the heater surface in the horizontal up orientation. In Fig. 23 with 
the surface in the vertical position, the heater surface temperature dropped to a quasi-steady 
boiling level after incipient boiling, which was considerably lower than that for the 
horizontal up orientation. This was expected for the heater surface in the vertical 
orientation, since the buoyancy force and the inertial force acted together to remove the 
vapor, which otherwise had an insulating effect, from the heater surfaces. As in the pool 
boiling vessel with the heater surface in the horizontal down position, a single vapor bubble 
covered the heater surface in the forced convection boiling loop with the heater surface 
horizontal down, even at the highest flow rate here, 4.5 cm/s. As a result, the heater 
surface temperature was unsteady and continually increasing, as shown in Fig. 24. It may 
also be noted that the measured surface temperature followed the one dimensional 
conduction prediction up to incipient boiling. Approximate adherence to the one 
dimensional conduction prediction is a consequence of 83 percent of the input energy going 
into the quartz substrate. The effect of orientation in earth gravity on the heater surface 
temperature at a q" T of 4 W/cm 2 as shown in Figures 22 through 24 is representative of the 

effect of the heater surface orientation on the measured heater surface temperature at the 
other heat flux levels. 

Figures 25 through 28 present the heater surface superheats at the onset of boiling as 
a function of q" T for the various orientations, subcoolings, and fluid velocities. Figures 25 

and 26 are for the low velocity, 1.7 cm/s. The heater surface superheats lie in a band 
between 10 and 40°C for surface Q-16 in Figure 25, positioned horizontal up, and between 
15 and 55°C for surface Q-13 in Figure 26, positioned horizontal down. As in the case of 
pool boiling, the time delay between the onset of boiling and the spread of boiling is most 
prevalent near a q" T of 8 W/cm 2 for surface Q-13. The largest heater surface superheats at 

the onset of boiling are for the low subcooling tests for a given heater surface orientation. 
The largest surface superheat at a velocity of 1.7 cm/s occurred for Q-13 in the vertical 
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position in Figure 26. In Figures 27 and 28 for the high velocity of 4.5 cm/s, the greatest 
heater surface superheats are for the low subcoolings in the horizontal down position. In 
Figure 28 for surface Q-16 at the fluid velocity of 4.5 cm/s, the heater surface superheats 
show no time delay between incipient boiling and boiling spread, and the superheats lie in a 
band between 25 and 40°C. From examination of Figures 25 through 28 for the vertical 
orientation, it may be concluded that the bulk velocity effects are small relative to buoyancy 
induced convection. 

Figures 29 and 30 present the delay time for the onset of boiling as a function of q" T 
with two subcooling levels with the horizontal down and vertical orientations for surface Q- 
13, for flow velocities of 4.5 cm/s and 1.7 cm/s, respectively. For the same subcooling 
level, the horizontal down surfaces tend to have smaller boiling inception delay times than 
the vertical heater surfaces, attributed to the assistance of buoyancy in the vertical 
orientation. The range in delay time before incipient boiling is nearly the same at both flow 
velocities. A time interval between incipient boiling and boiling spread, as was the case in 
pool boiling, occurred at the q" T of 8 W/cm 2 . As also was the case observed with pool 
boiling, the delay time to boiling inception decreased with increasing q" T , as shown in 
Figures 29 and 30. Figures 31 and 32 parallel Figures 29 and 30 for surface Q-16, but 
with the horizontal up in place of the horizontal down surface orientation. For the same 
subcooling, the horizontal up surfaces have larger boiling inception delay times than the 
vertical heater surfaces, and the range in delay time before the onset of boiling is essentially 
the same at both flow velocities. 

5.2 Quasi-Steady Boiling 

The early measurements of the forced convection boiling heat transfer behavior in the 
flow loop described here were made in two groups, designated A and B here, depending 
on the orientations of the heater surface relative to earth gravity, and are so presented 
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below. In all cases the liquid velocities are quite low, so that buoyancy indeed plays a role, 
and the liquid velocity vector is parallel to the heater surface as illustrated in Fig. 7. Thus, 
the angle between the flow velocity and gravity is changed as the loop is rotated, and is 
defined in Fig. 21. 

Group A consists of testing conducted with flow in the vertical up direction, so that 
the flow velocity enhances buoyancy with € = 180° in Fig. 21. Opposing sets of heaters, 
including both the metal heaters of Fig. 14 and the quartz-gold film heaters of Fig. 9, were 
installed in the test section of Fig. 6. Since these heaters were fabricated so as to be 
identical, operation in the horizontal orientation would simultaneously provide results for £ 
= 90° and B = 270° in Fig. 21, which are also included in Group A. However, as will be 
demonstrated, the behavior of the two surfaces presumed to be identical were not the same, 
so that results of changes between the vertical and the two horizontal orientations could 
only be compared in a qualitative manner. The differences in the changes were 
nevertheless quite dramatic. 

Subsequent to the testing which constituted Group A, the loop was modified so that it 
could be rotated through a total of almost 360°. In this way the same heater surface can 
now be subjected to all possible orientations between the flow direction and earth gravity. 
Results of these tests to date are presented as Group B below, and include only the use of 
the quartz-gold film heaters. 

In using the metal heater of Fig. 14 tire determination of the boiling heat flux from the 
measured power input requires that the heat loss from the lateral sides and the rear be 
estimated. Referring to the representative quantities plotted in Fig. 33, from the measured 
power input q^ with non-boiling convection on the fluid side, the computed heat transfer 
to the fluid, q FC , together with the computed heat transfer to the air q A using well- 
established correlations, the heat loss q^ is determined by the difference. Since the primary 
path of the loss cannot be identified precisely, an attempt was tnadg to relate this loss to the 
average of the measured front and rear surface temperatures, T w and T w b, respectively. 
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and to what is deemed to be the main sink for the heat losses, the test section housing 
whose temperature is close to the fluid temperature Tf. These are used to define the 
abscissa in Fig. 33. The ordinate on the left is expressed in terms of the power or heat 
transfer for the system of Fig. 14, while that on the right is expressed in terms of the heat 
flux at the boiling heat transfer surface area. It should be noted that the heat loss to the air 
out of the back side is negligible, and that the heat loss q L is virtually a linear function of 

the temperature difference defined above. Although the losses here appear to be a major 
part of the measured power input, the temperature differences are quite large, whereas with 
nucleate boiling taking place these differences will be considerably smaller, and the heat 
losses will be represented by the values on the left portion of the curves. Such calibration 
curves must be generated for each surface and for each orientation. 

5 2. 1 Group A Results: Metal and quartz heaters; 0 = 90°, 180°, 270° . 

Initial nucleate boiling results are shown in Figures 34 and 35 for the conditions 
given at the top of each figure. The numbers accompanying each data point represent the 
sequences in which the data were taken. The two sets of curves reproduce each other 
surprisingly well, and such should not normally be anticipated. Noteworthy are the large 
heater surface superheats prior to the onset of nucleate boiling, the lack of hysteresis effects 
in the nucleate boiling curve, and the distinctive cessation of boiling at heater surface 
superheats considerably lower than the onset. The non-boiling convection data at the end 
coincide with that at the beginning. The relatively large increase in heat flux once boiling 
begins occurs with a small increase in power, and is a result of the considerable decrease in 
heat loss accompanying the decrease in heater surface temperature. Subsequent to these 
tests the power supplies were modified to provide heat flux levels up to 20 W/cnA 

Fig. 36 shows the influence of orientation changes between vertical and horizontal-up 
and horizontal-down at the low velocity used. The non-boiling measurements for the three 
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orientations are indistinguishable from one another. In the vertical- upflow configuration, 
the two surfaces show a displacement of heater surface temperature by about 2 1/2°F at any 
given heat flux, to be contrasted with the data in Figs. 34 and 35, which were virtually 
identical. Rotation to horizontal-up (0 = 90°) produces an enhancement in performance, 
attributed to a decrease in the superheated boundary layer thickness because buoyancy is 
acting perpendicular to the flow direction. On the other hand, rotation to horizontal-down 
(0 = 270°) results in dryout of the surface, since buoyancy holds the vapor generated 
against the heater surface, and the velocity is too low to produce shear stresses sufficient to 
move the vapor away. The associated heat transfer rate is lower than the single phase 
liquid convection because of the lower thermal conductivity of the vapor, even with the 
reasonably high degree of inlet liquid subcooling. 

The influence of two low velocities (1.7 and 4.3 c/m s) are demonstrated in Fig. 37 
for the vertical orientation with liquid upflow, and in Fig. 38 for the horizontal orientation, 
facing both up and down. The change in heat flux with velocity in the non-boiling domain 
is clearly seen with the vertical up flow in Fig. 37, and with the horizontal orientation 
facing upward in Fig. 38. This is evidence that mixed convection, the combination of 
single phase free and forced convection is taking place. On the other hand, when the 
horizontal orientation is changed to facing down, also in Fig. 38, boiling and dryout again 
occurs, with no observable difference in behavior for the two velocities used. These two 
low velocities likewise have no influence on the nucleate boiling process in the horizontal 
orientation -facing upward (0 * 90°), which is indicative that buoyancy effects are 
dominating the process here. This is to be contrasted with nucleate boiling in the vertical 
upflow orientation in Fig. 37. For the low velocity of V * 1.7 c/ms the two heater 
surfaces differ by about almost 3°F in superheat, owing to somewhat different nucleating 
characteristics. When the flow velocity is increased to 4.3 cm/s, these now ta ke on an 
identical behavior, and with a larger slope. 
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Comparisons between the two types of heater surfaces used are presented in Fig. 39 
for the vertical orientation with upflow (0 = 180°) and in Fig. 40 for the horizontal 
orientation facing upward (9 = 90°) and downward (0 = 270°). 

For the vertical orientation in Fig. 39, differences between the opposing pairs of 
heaters are noted, with constant heater surface superheat difference of about 7°F for the 
quartz heaters and 2°F for the metal heaters. For this same velocity, the two metal heaters 
had identical surface superheats in the vertical orientation in Fig. 37, for which a liquid 
subcooling of 23.3°F was present, to be contrasted with 13.2°F in Fig. 39. This 
demonstrates that liquid subcooling can play a significant role. In spite of the larger 
differences between the quartz and metal heaters, the slope of heat flux-heater surface 
superheats are quite s imilar , with the metal surface superheats displaced to lower values, as 
is well known owing to its "better" nucleating characteristic. The negative slope of the 
quartz heaters at the lower levels of heat flux, with the data obtained with increasing heat 
flux, are manifestations of hysteresis effects, which are magnified considerably with 
smooth surfaces such as is present with a polished quartz surface. Stable nucleating sites 
are less likely to be present, and the spreading of the boiling process over the heater surface 
must be against both buoyancy and the flow direction. Hysteresis effects will be presented 
in somewhat more detail in the following section. As must be the case with the vertical 
orientation, the behavior of all four surfaces is identical in Fig. 39 in the non-boiling 
domain. 

Displacement of the heater surface superheat persists between the quartz and metal in 
the horizontal orientation facing upward (6 = 90°) in Fig. 40. A smaller hysteresis is 
evident with the quartz heater in this orientation, compared to the vertical one of Fig. 39, 
since the increase or spreading in the number of nucleation sites no longer needs to act 
against buoyancy, but only against the flow velocity. Only one data point was obtained 
with boiling with the quartz surface in the horizontal facing down orientation (0 = 270°), 
since the gold film used at this time could not withstand the relatively high surface 
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temperatures encountered. This point lies quite close to the data obtained with the metal 
heater, as it should since surface characteristics have no influence in the vapor bound or 
dry -out circumstance. 

5.2.2. Group B Results: Quartz heater; B = variable. 

The results presented below were obtained over a wider range of orientations than 
those in the previous section, with only the gold film quartz heater used. Since the 
spectrum of orientation angles covered is as yet incomplete, detailed interpretations of the 
results would be premature and will not be given here. 

Fig. 41 presents the complete data obtained with both increasing and decreasing heat 
flux, and it is noted that the behavior for these two cases is identical in the steady state with 
a sufficiently high level of heat flux, 6 W/cm 2 in this case. 

Figs. 42-46 demonstrate the influence of bulk liquid subcoolings of 4°F and 20°F 
for orientation angles varying ffom 0° to 180°, from flow directly counter to buoyancy to 
flow parallel to buoyancy. 

Data such as Figs. 42 - 46 are combined in Figs. 47 - 49 to illustrate the influences of 
orientation with the low velocity of about 4.3 cm/s used. Figs. 47 and 48 apply for low 
levels of subcooling, while Fig. 49 applies for the high subcooling. The increasing 
numbers on the plots correspond to increasing magnitudes of the orientation angle. All of 
the results in Figs. 42 - 49 were obtained with increasing heat flux. 

Fig. 47 includes only lower levels of heat flux, with orientation angles in the 
inverted" positions, while Fig. 48 has higher levels of heat flux, with orientation angles in 
the upward-facing positions. 

Fig. 49 covers a wider range of orientations and levels of heat flux with the high 
subcooling of 20°F. 
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L FU TU RE WORK 

a. Examination of the steady state nucleate boiling results with low velocities, in 
section 5.2 above, makes it clear that in conducting such experiments in earth gravity little 
can be discovered about the potential behavior for corresponding conditions in long term 
microgravity. Experimental work in both short and long-term microgravity will be required 
to understand the basic elements that constitute forced convection nucleate boiling at the 
smaller flow rates with low Reynolds Numbers, and should include heat flux levels up to 
the critical values. 

b. A number of extensionsof the initial work presented here can be carried out in 
earth gravity: 

i. Operation at near saturated liquid at the entrance to the test section, in which 
it is anticipated that the influences of orientation will be more dramatic. 

ii. Extend the levels of heat flux up to the critical or "burn-out" level as 
functions of orientation, subcooling, and low-level fluid velocities. 

iii. Study the influence of small variations from the horizontal down orientation 
(0 = 270°) with various small velocities to establish the conditions necessary to remove the 
vapor from the heater surface. 
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Fig. 1. Influence of Buoyancy on Temperature Distribution 
in Vicinity of Heater Surface with Pool Boiling. 
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Lift 



Fig. 2. Forces Acting on a Growing/Collapsing Vapor Bubble Attached to a Wall. 
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Thermal Management in Space 

(Change in pumping power requirements using phase change) 

rh 

Qout 

( Radiation 
to space) 


liquid heat capacity : 
latent heat : 

or : 

Power : AP « / • V 2 • L in laminar flow 



Also: Mass of fluid required is reduced by 1/2 !!! 


Fig. 3. Example of Advantage in use of Phase Change to Transport Energy. 
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Fig. 4. Early Flow Loop for the Study of Forced Convection 

Boiling of Various Orientations and at Reduced Gravity. 
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Note : Ml. M2 & M3 — Cold Mirrors 
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Fig. 8. Optical Paths for Two View Photography. 
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Fig. 9. Gold Film Heating Surface as Simultaneous Resistance Thermometer. 
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TEST FLUID 


Tsat 



- TRANSPARENT GOLD FILM 
(HEATER + RESISTANCE 
THERMOMETER) 


Fig. 10. Demonstration of Concept for Transient Heating with a Thin Film. 
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Time (Sec) 


Fig. 1 1. Thin Gold Film Temperature with Onset of Heating in Pool Boiling at a/g = 1 




Fig. 12. Thin Gold Him Temperature with Onset of Heating in 
Pool Boiling at a/g = 10*^ for 5 seconds. 

49 




Fig. 13. Heat Transfer Between Two Communicating Semi-infinite 
Solids with Uniform Plane Heat Source at the Interface. 
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Fig. 14. Copper Heater Surface for Quasi-steady Operation. 


51 




Fig. 15. Pressure Control Bellows Assembly. 
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HEISE GAUGE MODEL 623 (or Socrm-Hodol Z70) 
0- SOPS! A 
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Fig. 16. Flow Loop Pressure Level Control System. 


53 









L. 

Q) 



54 


Rg. 17. Schematic of Loop Row Control System. 







Mein Loop Pressurizer 


R 1 1 3 

3 Pass 
Counterflov 
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Fig. 18. Schematic of Loop Condenser/Cooler Control System. 





Fig. 19. Schematic of R-l 13 Degassing Unit. 
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Nucleation & 
Incipient Boiling 


Forced-Convection 
Nucleate Boiling 



Transient 



Quasi-Steady 


Pool 


Forced Convection 


Steady-State 

Boiling 


Independent Variables 

Fixed Fluid R-113 

Fixed Pressure (near Ambient) 


Heater Surfaces 


Semi-transparent Gold on Quartz 
^Gold on Metal (Copper) 


Variable AT subcooling 

Variable Velocity (velocity distribution) 
Variable Heat Flux 
Variable Orientation 


Dependent Variables 

Conditions for incipient boiling 

Vapor bubble motion 

Mean heat transfer behavior 


Fig. 20. Experimental Parameters for Flow Loop Operation in Earth Gravity. 
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Inlet 

Flow 


Orientations used in experiments 



Fig. 21. Flow Gravity Boiling Heat Transfer with Various Earth Gravity Orientations. 
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Fig. 22. Thin Film Gold Heater Temperature in Forced Convection 
Boiling with Horizontal Up Orientation. 
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Fig. 23. Thin Film Gold Heater Temperature in Forced 
Convection Boiling with Vertical Orientation 
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Fig. 24. Thin Him Gold Heater Temperature in Forced Convection 
Boiling with Horizontal Down Orientation. 
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Fig. 25. Heater Surface Superheat at Incipient Boiling for Surface Q- 16 at 1.7 cm/s. 
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Total Heat Flux, q*r (W/cm 2 ) 


Fig. 26. Heater Surface Superheat at Incipient Boiling for Surface Q-13 at 1.7 cm/s. 
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Fig. 27. Heater Surface Superheat at Incipient Boiling for Surface Q-13 at 4.5 cm/s. 
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Fig. 29. Incipient Boiling Delay Time at 4.5 cm/s for Surface Q-13. 
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Total Heat Flux, q"x. (W/cm 2 ) 


Fig. 30. Incipient Boiling Delay Time at 1.7 cm/s for Surface Q-13. 
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Fig. 32. Incipient Boiling Delay Time at 1.7 cm/s for Surface Q-16. 
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40 


Data based on Tf = 107. 3°F 



Fig. 33. Representative Heat LossCalibration - Metal Heater No. 2 - Vertical. 
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q/A - w/cm 



q”(w/cm2) 


P»l7.91p«ii. Tmi»12*.6*F, AT*u^ 18»F. r»1.72cm/i 



Fig. 34. Sample Forced Convection Nucleate Boiling Results. Subcooled, 
Vertical, Low Velocity, Metal Heater No. 1. Test No. 62289. 
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(M7.91ptu, Tmib 12S.6*F. AT»u)^18*F, r»1.72an/i 



Fig. 35. 



ATwsup(K) 


Sample Forced Convection Nucleate Boiling Results. Subcooled, 
Vertical, Low Velocity. Metal Heater No. 2. Test No. 62289. 
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Fig. 36. Influence of Orientation with Forced Convection Nucleate Boiling. 
Metal Heaters. ATsub = 23.2°F, V = 1.7 cm/s. 
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Fig. 37. Influence of Velocity with Forced Convection Nucleate Boiling and Vertical 
Up How 180°). Metal Heaters. AT** » 23.3°F. 
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Fig. 38. Influence of velocity with forced convection nucleate boiling and horizontal 
orientation. Up (0 = 90°) and down (&•= 270°). Metal heaters. 

ATsub = 13.FF. 
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Tw*up (K) 


Fig. 39. Comparison of metal heater versus gold film quartz heater with otherwise 
identical conditions: vertical up flow V ($= 180°); ATsub = 13.2°F; 

P - 16 psia; V = 4.3 cm/s. 





q"(W/cm2) 



Fig. 40. Comparison of metal heater versus gold film quatrz heater with otherwise 
identical conditions: Horizontal facing up (0= 90°) and down (0-= 270°); 
ATsub = 13.1°F; P - 16 psia; V = 4.3 cm/s. 
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Fluid: R-113 
Surface Q23 
Rowrita « 4 . j e 'y / S 
Subcooiing ■ 5 F 
Inlet temp ■ 1 17 F 
- 0 - - 225 degrees 


Q Increasing flux 
* Decreasing flux 


Fig. 41. Hysteresis Effects with quartz heater at = 225°. Test No. 112090. 



Heat Flux (W/cm*2) 


Heat Flux vs. Superheat 



Angle - 0 degrees 
Inlet Temperature « 120.0 
Flow Rate ■ 4.3 cm/s 
Gold Heater Surface (Q23) 


n ATsub « 4 F 
• ATsub ■ 20 F 


AT superheat (F) 

Fig. 42. Influence of Subcooling for-0- = 0°. Test No. 022591. 
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Heat Flux (W/cm A 2) 


Heat Flux vs. Superheat 



Angle « 90 degrees 
Inlet Temperature * 120.0 
Flow Rate ■ 4.3 cm/s 
Gold Heater Surface (023) 


D ATsub ■ 4 F 
• ATsub = 20 F 


AT superheat (F) 

Fig. 44. Influence of Subcooling for# = 90°. Test No. 022591. 
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Heal Flux (W/cm^2) 


Heat Flux vs.' Superheat 



Angle ■ 135 degrees 
Inlet Temperature = 120. 
Flow Rate ■ 4.3 cm/s 
Gold Heater Surface (02: 


o ATsub m 4 F 
♦ ATsub = 20F 


AT superheat (F) 

Fig. 45. Influence of Subcooling for O' = 135°. Test No. 022591. 
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Heat Flux vs. Superheat 



Angle ■ 180 degrees 
Inlet Temperature - 120.0 F 
Flow Rate ■ 4.3 cm/s 
Gold Heater Surface (Q23) 


B ATsub ■ 4 F 
e ATsub « 20 F 


AT superheat (F) 

Fig. 46. Influence of Subcooling for -O- = 180°. Test No. 022591. 
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Fig. 47. Influence of Orientation with Low Subcooling. 
Low Heat Flux. Test No. 112090. 
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Heat Flux (W/cm A 2) 


Heat Flux vs. Superheat at Varying Orientations 



Inlet Temperature = 120.0 
Subcooling * 4.0 F 
Flow Rate ■ 4.3 cm/s 
Gold Heater Surface (Q23) 


AT superheat (F) 


□ Angle = 


• Angle « 4 
( 3 ) ■ Angle -9 

♦ Angle = 1 
■ Angle • 1 


Fig. 48. Influence of Orientation with Low Subcooling. Test No. 022591. 
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Appendix A 


Finite Difference Solution to 
Mixed Convection in a Rectangular Duct 


by 

Kevin M. Kirk 
December 19, 1990 


A two-dimensional transient finite difference solution is presented 
for the problem of laminar mixed convection in a rectangular duct. The 
incoming flow to the duct, which can be specified as being either uniform 
flow or fully developed poiseuille flow, passes over a heated patch 
subjected to a constant heat flux. The problem is solved for varying 
orientations but with particular attention given to vertical upflow, vertical 
downflow, and horizontal flow. 
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I. Introduction 


The presence of gravity can give rise to significant buoyancy forces in low 
speed flows where behavior predicted solely from the consideration of 
forced convection can lead to erroneous results. Furthermore, the direction 
of gravity relative to the heating surface can have significant influence on 
the importance of gravity in natural convection. 

Nyugen et al. (1979) studied mixed convection for the case of horizontal 
flow between parallel plates with isothermal walls and nonuniform 
velocity and temperature profiles. Their results demonstrated that heat 
transfer above a heated place was enhanced while that below a heated 
plate decreased with respect to purely forced flows. Other analytical 
investigations studying the effect of buoyancy on fully developed laminar 
flow between parallel planes (Gill and Del Casal, 1962) have showed that a 
separated recirculating flow can occur at the upper wall. Kennedy and 
Zebib (1982) presented numerical and experimental work on the effects of 
free convection on forced laminar flow with a local heat source on the 
lower wall. Their results demonstrated that a recirculation region adjacent 
to the top wall and above the heat source was generated due to the effects 
of natural convection. A numerical study by Davis (1976) examining the 
effect of orientation on mixed convection over a heated patch of constant 
temperature predicted recirculation for the case of vertical downflow. 

The present numerical study was done with the assumption of constant 
heat flux on a small heated patch. The heat transfer characteristics and 
flow patterns were studied for four different orientations where particular 
attention is given to the occurrence of flow recirculation since it is 
well-known that an inflection in velocity can give rise to instabilities. 

A study of mixed convection also has relevance in the study of incipient 
forced convection boiling for the case of low speed flows. For example, 
when a sufficiently high heat flux is applied to a surface, boiling will 
commence when the surface has locally attained the required superheat 
for boiling inception. Knowledge of the temperature distribution in the 
vicinity of the heater surface is critical to understanding the initial 
nucleation and subsequent spreading phenomenon. Hence, if buoyancy is 
assumed significant, a transient mixed convection solution is necessary and 
such a solution is given numerically in this work for laminar flow. 
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II. Problem Formulation 


A. Dimensional Governing Equations 

The geometry of the problem is illustrated in Figure 1 with the inlet flow 
being directed in the +x direction and horizontal flow with the surface 
facing upwards occuring when 0 = 90°. 



Figure 1 - Two dimensional problem geometry 


The standard problem of two-dimensional mixed convection in a 
rectangular field containing a Newtonian fluid for which the Boussinesq 
approximation may be applied can be described by the following 
dimensional equations: 


Continuity 


du dv 
dx dy 


0 


( 1 ) 


Momentum 

in the x-direction 


du.du.3u i 3P 

-=r- + Ux- + V^r- = — l -5~ 

dt dx dy p 0 dx 


-gpCr-Tjcose+i 


ft ft' 

dx 2 dy 2 . 


( 2 ) 
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M omentum 

in the y-direction 


dv , dv , 9v i dP 

dt dx By p 0 dy 


+ gp(T - T 0 ) sin0 + 


Po dx 2 By 2 


Energy 


dT dT dT 

dt dx dy 


a 


&T ( d*T 
dx 2 By 2 


(3) 


(4) 


Alternatively, the problem can be formulated in terms of the stream 
function and vorticity. This approach ensures automatic satisfaction of the 
continuity equation, eliminates the pressure gradient term and reduces the 
number of differential equations to be solved by one. 


The stream function, V, and vorticity, to, are defined as 


dw dw 

V "’d7 


dv du d^V dfy 

dx dy dx 2 3y2 


(5) 


Differentiating equation (2) with respect to y and equation (3) with respect 
to x and, then, subtracting (3) from (2) yields the following equation after 
making appropriate substitutions for vorticity and the stream function 


Bat dydoo dydco 
dt dy dx dx By 



dT 

dy 


:os6 


LiL 

d 2 © 

rpo 

dx 2 



( 6 ) 


After substituting the stream function into the energy equation for the 
velocity, the energy equation becomes 

dT dydT dydT 
dt dy dx dx dy 


d 2 T | B 2 T 
dx 2 dy 2 , 


Equations (5), (6), and (7) form a set of three coupled non-linear 
differential equations to be solved numerically with dependent variables 
T, a), and V. 
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B. Boundary and Initial Conditions 

From the governing equations, it is observed that one condition for each of 
the variables, T, to, and V, are needed on each boundary as well as an initial 
condition for each variable. 


Boundary conditions on V at y=0 and y=L are obtained from the no slip and 
no penetration conditions. That is, on a solid boundary, V = constant, and 
= 0, where *1 is the direction normal to the solid surface. The 
boundary condition for V at the inlet (i.e. x=0) is obtained from specifying 
either poiseuille or uniform flow as the inlet flow condition. At the exit, it 
is assumed that the curvature of the streamlines is negligible. Hence, at 
the exit plane 



x«xii,y,t) — o 


Boundary conditions for vorticity, to, must be expressed in terms of the 
stream function. For uniform flow at the inlet, x=0, vorticity is allowed to 
diffuse upstream. Thus, applying equation (5) at the inlet we get 


0 ) = - 


dty 


For the inlet condition of poiseuille flow, the boundary condition at x=0 is 
obtained simply by applying equation (5) at the boundary. 


Along the wall, y=0, the flow is assumed to be essential parallel, hence, 
dy/3x = 0 and 3/Sx = 0 From these assumptions, the following second-order 

Woods formula can be obtained through application of equations (5) and 
(6) at an insulated wall for steady flow 


o>u= + 


3(yu- Vi.2) 
(Ay) 2 


Finally, at the exit plane since the curvature of the streamlines is 
considered negligible, equation (5) reduces to 


<o = * 


dy2 
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The boundary conditions for temperature are a constant inlet temperature, 
adiabatic wall surfaces, a constant heat flux applied at the heater surface, 
and negligible curvature of the isotherms at the exit. 


The initial and boundary conditions are summarized as follows: 
at t=0 

T(x,y,0) = T„ 
for uniform flow 

y(x,y,0) = U„y 
co(x,y,0) « 0 

for poiseuille flow 

V(x,y,0) = 6Uj^. 

du 


3L 2 J 


<o(x,y,0) = ~ = -6U„ 


X.3J£] 

L L 2 


at x=0 


T(0,y,t) = T„ 


for poiseuille flow 

V(0.y,,).«U^-^ 

«* 

for uniform flow 
V(0,y,t)=U - y 

®(0,y,t) = 

ox 2 


at x=x ex jt 


-r-HXexit,y,t) * 0 
ox 2 
j2t 

- 0 

ox 2 


w(x«xii.y»t) 


9V 

dy 2 
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at y=0 


V(x,0,t) = 0 

fo'.o.t) = 0 

O)(x,0,t) = - ;ko 

2 i.2 


, 3(Vu- VU2) 
(Ay) 2 


on heater surface 

f^x,0.t) = -£ 
dy k f 

at y=L 

3 “ (x,L,t) = 0 
dy 

V(xO-,t)=U«L 

0Xx.L. t ) - • h> 

2 i-M’ 1 (Ay) 

C. Non-Dimensionalization of Equations 

It is usual to put the equations into dimensionless form to reduce the 
number of independent parameters and, hence, the computational effort. 
Non-dimensionalization also offers opportunities for enhanced physical 
understanding through the representation of the relative forces into 
dimensionless parameters. 


In order to make the equations (5) thru (7) dimensionless, the following 
dimensionless variables are defined where the characteristic velocity is the 
mean flow velocity, U«, and the characteristic length is the heater length, 1. 


5T — x/1 y * y/1 


Y=V/(UJ) 


o)«co(lAJJ 
~ T-T„ 

e “(qTO 


t 


mtSL 

l 2 


Substituting the dimensionless variables into the governing equations, the 
following dimensionless equations are obtained 


I dco dydS) dydu ) b 
P e + dy 95c dx dy 


D . 90 0 d8. 0 

Ri Ur^cosO + T=sin0 
\dy dx J 



( 8 ) 
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( 9 ) 


( 10 ) 


illustrated in Figure 2. 




For polseullle flow: 


VOUAO-L/l 


L ,rj c (AJO 




For uniform flow: 

St°.y,t) - ■ ^ 

it 1 


t L 
•» 


EJ 2 


a ~ . aKr 

3 _ ay* 

-> a 2 * 

-» ~<i^,y,t)-o 
■* ex 4 


Y<o,yt) - y 




<A 30 


V(r,0,t)-0 

a) Vorticity and stream function boundary conditions 




e<o,y,t)-o ^ l 


30 

^Or,L/i,t)-o 


t ^|<5r„*,y,t)-0 
■» axr 2 


as A 

ay*™- 0 fer.o.'«-i 

b) Temperature boundary conditions 


Figure 2- Dimensionless Boundary Conditions 
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III. Method of Solution 


The second order derivatives were approximated using central differences. 
However, approximating the convection terms in the governing equations 
with central differences suffers from the well known Reynold’s cell 
number restriction. This instability may be overcome with the use of 
upwind differencing, instead of central differencing, since diagonal 
dominance is ensured. However, it has been recognized (Davis, 1974) that 
the first order upwind difference approximation can generate significant 
truncation errors which may give rise to a "false diffusion," especially in 
the case of high Reynolds number flow where physical diffusion is weak. 

Nonetheless to assure convergence over a wide range of Reynolds number 
an upwind differencing scheme was used where the local stream direction 
was determined from the the stream function values around the point in 
question. 

The equations were solved using an alternating direction implicit method. 
The method employs two difference equations which are used in turn over 
successive time steps, each of duration At/2. The first equation is implicit 
only in the x-direction, while the second equation is implicit only in the 
y-direction. The advantage of this method is that this technique requires 
only the solution of a tridiagonal matrix which affords a straightforward 
solution. 

A stability analysis was performed experimentally through modification of 
the mesh size and time step until convergence was achieved for a given set 
of dimensionless parameters. 
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IV. Results and Discussion 

Representative contour plots are given of the isotherms and the stream 
function in Figures 3 thru 9 for different orientations and values of Ri, and 
Re. In all the figures it is assumed that Pr=l, hence. Re = Pe. Low values 
of Re were used in the representative contour plots in order to accentuate 
the effects of buoyancy. 

It is seen from the contour plots of vertical upflow in Figure 6 that the 
buoyancy acts to accelerate the flow near the heated surface thereby 
enhancing heat transfer. The acceleration of the flow is illustrated by the 
increased density of the streamlines near the heating surface, and is 
accompanied by an increase in the mean Nusselt number across the heater 
as shown in Figure 10. In the contour plot of vertical downflow, shown in 
Figure 4, the flow is decelerated due to the buoyant force acting up on the 
heated fluid near the heater surface, and the mean Nusselt number across 
the surface decreases. 

If the Richardson number is sufficiently large in vertical downflow, the 
numerical solution predicts flow reversal and circulation where the point 
at which the onset of circulation occurs is shown in Figure 10. 


Inlet Poiseuille Flow 
Vortical Orientation 
Pr- 1 


■ Re-1 0 Vertical Down 
e Re-30 Vertical Down 

■ Re-100 Vertical Down 
e Re-10 Vertical Up 

■ Re-30 Vertical Up 
o Re-100 Vertical Up 

Onset of Circulation 


0 10 20 30 40 

Ri 



Figure 10 - Mean Nu vs Ri for Vertical Flow 
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Without circulation, the steady-state maximum temperature in vertical 
downflow occurs at the downstream edge of the heater surface, however, 
when circulation does occur, the steady-state maximum temperature is 
located at the upstream edge of the heater surface. 

It is interesting to note that for a range of Re from 0 to 100, it is shown 
that the value of Ri at which the onset of circulation occurs seems to be 
solely a function of Re for a given Pr, and is defined below as the critical 
Richardson number, Ri cr . In other words, as Re increases, less buoyancy is 
required for the onset of recirculation. This trend has also been predicted 
by Davis (1976). This relationship is shown graphically in Figure 11, and is 
given by the equation 

Ri„ = 48.98Re* 0 368 

or in terms of the Grashof number 

Gr cr = 48.98Re 1632 



Vertical Downflow 


Inlat Poisauilla Flow 
Pr - 1 


■ Critical Ri 


Figure 11 - Critical Ri vs Re for Vertical Flow 

In horizontal flow with the surface facing up, the streamlines initially bend 
toward the heated surface as shown in Figure 8 indicating that the flow 
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upstream of the surface is drawn toward the surface. However, at near the 
mid-point of the heater surface the flow begins to move away from the 
surface due to the influence of buoyancy. Conversely, in horizontal flow 
with the surface facing down, illustrated in Figure 10, the streamlines 
initially bend away from the surface and, again, at near the mid -point of 
the heater the flow is drawn back toward the heater side due to buoyancy 
effects. 

As expected, the highest temperature on the heater surface in horizontal 
flow was always at the downstream edge of the surface as seen in the 
contour plots of the isotherms, Figures 7 and 9. Furthermore, a higher 
maximum temperature was attained with the surface facing down than 
with the surface facing up in horizontal flow with identical values of Ri, Re, 
and Pe. The higher temperature is due to the degradation of the heat 
transfer with the surface facing down as illustrated by the decrease in the 
Nusselt number with increasing Richardson number in Figure 12. 

On the other hand, buoyancy enhanced heat transfer when the surface was 
facing up in horizontal flow which is also demonstrated in Figure 12. 

Recirculating flow is also predicted for the case of horizontal flow with the 
surface facing downward for large values of Ri, however, whether or not 
this represents reality remains a question left to experimentation. 



Inlet Poiseuille Flow 
Horizontal Orientation 
Pr - 1 


■ Re-30 Horiz. Up 
e Re-50 Horiz. Up 

■ Re-100 Horiz. Up 

♦ Re-30 Horiz. Down 

■ Re-50 Horiz. Down 
a Re-100 Horiz. Down 

^ Onset of Circulation 


Figure 12 - Mean Nu vs Ri for Horizontal Flow 
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The effect of buoyancy on heat transfer at orientations other than 
horizontal and vertical orientations is given in Figure 13 for a constant 
Reynold's number. 


R« - 1 00 
Pe - 100 

1 ■ Nu-0 

2 ♦ Nu-45 

3 ■ NU-90 

4 ♦ Nu-135 
SB Nu-180 

6 ° Nu-225 

7 * Nu-270 

8 a Nu-80and238 


0 1 2 3 4 5 6 7 8 9 10 1 1 12 13 14 IS 16 

R1 

Figure 13 - Mean Nu vs. Ri at Varying Orientations 

As 6 increases from 0° (vertical downflow), there is a considerable 
decrease in the tendency for the buoyant force to degrade the heat 
transfer. This trend continues until G equals $ 1 , where Ri has little 
influence on Nu and can be approximated by a horizontal line as seen 
on Fig. 13. By interpolation, the angle of $1 has been determined as 
approximately 80 degrees. As 6 increases past $1 the buoyant forces 
no longer degrade the heat transfer from the plate but instead help 
to enhance it by increasing Nu. The greatest value for Nu occurs 
between 6 = 135 and 6 = 180 since components of buoyant forces and 
forced convection are in the same directions for upward flow. 

As 0 rotates from 180 to 270 degrees there once again begins a 
degradation of heat transfer as the buoyant forces become more directed 
toward the heating surface. This ultimately results in warm fluid 
remaining closer to the heating surface and a decrease in heat transfer. 
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Figure 3 - Dimensionless Isotherms 
Vertical Downflow 
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Figure 4 - Dimensionless Scream Function 
Vertical Downflow 
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Figure 5 - Dimensionless Isotherms Vertical Upflow 
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Figure 6 - Dimensionless Stream Function 
Vertical Upflow 
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Figure 7 - Dimensionless Isotherms 

Horizontal Flow with Surface Facing Up 
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R 


100 


Re « 30 


Pe * 30 



Figure 8 - Dimensionless Stream Function 

Horizontal Flow with Surface Facing Up 
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Figure 9 - Dimensionless Isotherms 

Horizontal Flow with Surface Facing Down 
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V. Program Listing 


C PROGRAM TRANSADI 
C 

C BY KEVIN M KIRK THE GREAT 
C 

C This program solves for the transient temperature and velocities 

C fields in the presence of combined natural and forced 

C convection in a rectangular duct. The problem is solved 

C for the case of either slug or poiseuille flow over a heated patch 

C subjected to a constant heat flux. 

C 

C The program has the capability to solve 
C the problem for varying orientations. That is, gravity can be- 

C perpendicular or parallel to the heating surface. 

C 

C DECLARATION OF VARIABLES 

REAL*8 PSI( 1 50, 1 50), W( 1 50, 1 50),TH( 1 50, 150),A3»C,DX 
REAL*8 Q,UINF,G,TC,DIFFU,MU,RI,REJ’E4 J J > R,UP,UM,VP,VM 
REAL*8 E 1 ( 1 0000),PSIOLD( 1 50, 150),DIFF1 ,DIFF2,TINF,NU,NUT 
REAL*8 F,DT,A1(10000),B1(10000),C1(10000),D1(10000),NUX 
REAL*8 TOLD( 150,1 50), ALPHA^I 
INTEGER N^ 1 ^ 2 MU,SEIECT,nER,ORIENT,COUNT^STEP 
INTEGER NUMJC 
LOGICAL CONVERG 
C 

C PSI=DIMEN SIONLESS STREAM FUNCTION 

C W = DIMENSIONLESS VORTICITY 

C TH = DIMENSIONLESS TEMPERATURE 

C DX = MESH SIZE 

C TINF = TEMPERATURE OF FLUID AT INLET 
C Q = IMPOSED HEAT FLUX ON SURFACE 
C UNIF = INLET FLUID VELOCITY 

C G = GRAVITY 

C TC = THERMAL CONDUCTIVITY 

C DIFFU = THERMAL DIFFUSIVITY 

C MU = KINEMATIC VISCOSITY 

C RI = RICHARDSON NUMBER 

C RE = REYNOLD S NUMBER 

C PE = PECLET NUMBER 

C N = NUMBER OF NODES IN X DIRECTION 
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C M = NUMBER OF NODES IN Y DIRECTION 
C A = LENGTH OF CHANNEL 

C B = HEIGHTH OF CHANNEL 

C C = ENTRY LENGTH 

C ALPHA=ANGLE OF INCLINATION WITH VERTICAL 
C 

C Define default properties for freon (R-113) 

C 

DIFFU=4.5E-6 
MU=4.63E-7 
TC=0.066 
PI=3. 14159 

Initialize counters and time. Open output files. 

NCOUNT=0 
TIME=0.0 

OPEN(UNIT=4,FTLE='thcta.dat',STATUS=’NEW) 
OPEN(UNlT=5,FTLE='stream.dat , ,STATUS='NEW’) 


Input dimensions of duct and time step 


WRITE(*,200) 
READ(*,*)A,B.C 
WRITE(*,214) 
READ(*,*)L,DX 
WRITE(*,211) 
READ(*,*)DT,STOPJ > STEP 
NPRINT=NINT(PSTEP/DT) 
NSTEP=NINT(0.2/DX) 
N=NINT(A/DX)+1 
M=NINT(B/DX)+ 1 
WRITE(4,300) 
WRITE(5,295) 
N1=NINT(C/DX)+1 
N2=NINT((C+L)/DX)+ 1 
NUM=(N-2)*(M-2) 
I2=(Nl+N2)/2 

CONVERT TO METERS 


A= A/1 00.0 
B=B/100.0 
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C=C/100.0 
L=L/ 100.0 
DX=DX/100.0 
WRITE(*,202) 

WRITE(*,203) 

WRITE(*,204) 

READ(V)SELECT 

WRITE(*,208) 

READ(*»*) ALPHA 
WRITE(*,209) 

READ(*,*)NTYPE 

WRITE(*,212) 

WRITE(*,213) 

READ(*,*)MODE 
IF (SELECT.EQ.2) THEN 
WRITE(*,210) 

READ(*,*)TINF,BETA,Q 

WRITE(*,220) 

READ(V)UINF,G 

WRITE(*,230) 

READ(*,*)TCJ)IFFU 

WRITE(*,240) 

READ(V)MU 

MAKE UNITS CONSISTENT IN SI UNITS 

QssQ*(100.0**2) 

UINF=UINF/ 100.0 
G=G*9.8 

RI=(G*Q*(L**2))*BETA 

TEMP=((UINF**2)*TC) 

RI=RI/TEMP 

RE=UINF*L/MU 

PR=MU/DEFFU 

PE=RE*PR 

WRITE(*^15)RI,RE,PE 

ELSE 

WRITE(*,250) 

READ(*,*)RI,RE t PE 

ENDIF 

WRITE(*,256)N t NU42,M t NSTEP 
SET UP FILE FOR PLOTTING PROGRAM 


125 




c 


WRITE(4,257)((N-1 )/NSTEP)+l ,((M- 1 )/NSTEP)+l 
WRITE(5,257)((N- 1 )/NSTEP)+l ,((M- 1 )/NSTEP)+l 
WRITE(4,259)(((I-1 )*DX* 100),I= 1 .N^fSTEP) 

WRITE(5,259)(((I- 1 )*DX* 100)4=1 ,N,NSTEP) 

WRITE(4,259X((I- 1 )*DX* 1 00),I= 1 ,M,NSTEP) 

WRITE(5,259X((I- 1 )*DX* 1 00),I= 1 ,M ,NSTEP) 

PAUSE 

DX=DX/L 

C 

C SET INITIAL CONDITIONS FOR ALL NODES 
C ASSUME POISEUILLE OR SLUG FLOW, AND UNIFORM INITIAL 
C TEMPERATURE 
C 

DO 10,I=1,N 
DO 204=1, M 
Y=(J-1)*DX 

IF (NTYPE.EQ.l) THEN 
PSI(IJ)=Y 
W(IJ)=0.0 
ELSE 
S=L/B 

PSI(I,J)=6.0*((S*(Y**2)/2.0)-((S**2)*(Y**3)/3.0)) 

W(IJ)=-6.0*(S-(2.0*Y*S**2)) 

ENDEF 
TH(I4)=0.0 
20 CONTINUE 
10 CONTINUE 
C 

C INCREMENT TIME AND COUNTER FOR PRINTING OUTPUT 
C TO SCREEN 
C 

25 TIME=TIME+DT 
NCOUNT=NCOUNT+ 1 
C 

C UPDATE BOUNDARY CONDITIONS 
C 

C AT INLET ASSUME UNIFORM OR POISEUILLE FLOW, AND UNIFORM 

TEMPERATURE 

C 

1=1 

IF (NTYPE EQ.l) THEN 
DO 304=1, M 
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TEMP=2*PSI(I,J)-5*PSI(I+1 J)+4*PSI(I+2,J)-PSI(I+3J) 
W(I,J)=(-1.0/(DX**2))*TEMP 

30 CONTINUE 
ELSE 

DO 3U=2,M-1 

TEMP 1 =2*PSI(U)-5 *PSI(I+ 1 J)+4*PSI(I+2J)-PSI(I+3,J) 
TEMP2=PSI(IJ+ 1 )+PSI(U- 1 )-2*PSI(U) 

W(IJ)=(-1 .0/(DX**2))*(TEMPl +TEMP2) 

31 CONTINUE 

TEMP=2*PSI(I,1)-5*PSI(I,2)+4*PSI(L3)-PSI(I,4) 

W(l,l)=(-1 .0/(DX**2))*TEMP 
TEMP=2*PSI(I,M)-5*PSI(I,M- 1 )+4*PSI(I,M-2)-PSI(I,M-3) 
W(I,M)=(-1.0/(DX**2))*TEMP 
END IF 
C 

C AT Y=0, assume adiabatic surfaces everywhere 

C except on heater surface. Use Woods approximation for 

C vorticity on the surface. 

C 

J=1 

DO 40,I=2,N1-1 

TH(IJ)=( 1 .0/3 .0)*(4*TH(IJ+ 1 )-TH(I J+2)) 
W(U)=-0.5*W(U+1)-(3*(PSI(U+1)-PSI(I,J»/(DX**2)) 

40 CONTINUE 
DO 50,I=N1,N2 

TH(I J)=( 1 .0/3 .0)*(4*TH(I J+ 1 ) : TH(I J+2)+2.0*DX) 
W(U)=-0.5*W(I,J+1H3*(PSI(U+1)-PSI(U))/(DX**2)) 

50 CONTINUE 

DO 60,I=N2+1,N 

TH(IJ)=(1.0/3.0)*(4*TH(U+l)-TH(U+2)) 

W(U)=-0.5*W(U+1)-(3*(PSI(I,J+1)-PSI(U»/(DX**2)) 

60 CONTINUE 

C 

C At exit, assume constant temperature gradient, constant 

C curvature of streamlines 
C 

I=N 

DO 69,J=2,M-1 

TH(I J)=0.5 *(5 *TH(I- 1 r J)-4*TH(I-2 J)+TH(I-3 J)) 
W(UM-1.0/(DX**2))*(PSI(U-l)-2*PSI(U)+PSiaJ+l)) 
69 CONTINUE 

C 

C AT Y=B 
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c 

J=M 

DO 70,I=2,N 

TH(I,J)=(1.0/3.0)*(4*TH(U-l)-TH(U-2)) 

W(IJ)=-0.5 * W(I,J- 1 )-(3 *(PSI(U- 1 )-PSI(U))/(DX**2)) 

70 CONTINUE 
C 

C PRINT OUT TEMPERATURE AND STREAM FUNCTION TO THE SCREEN 
C AT PRESELECTED TIME INTERVALS 

C 

IF ((NCOUNT.EQ.NPRINT).OR.(TIMEUE.DT)) THEN 
WRITE(*,261)TIME 
DO 32J=M,1,-NSTEP 

WRITE(*,270)(PSI(I4)4=1 ,N .NSTEP) 

32 CONTINUE 

WRITE(*,201) 

DO 34J=M,1,-NSTEP 

WRITE(*,280)(TH(U) J=1 .N^STEP) 

34 CONTINUE 

WRITE(*,290) 

READ(*,*)KEV 
IF (KEV.EQ.1) THEN 
GOTO 125 
ENDIF 
NCOUNT=0 
ENDIF 
C 
C 

C USE ALTERNATING IMPLICIT DIRECTION SCHEME TO SOLVE FOR 
TEMPERATURES 

C AT INTERIOR NODES. START WITH THE X-DIRECTION BEING IMPLICIT 
C AND USE UPWINDING SCHEME BASED ON LOCAL VELOCITY DIRECTION TO 

C AID CONVERGENCE 
C 

DO 804=2,M-1 
DO 854=24^-1 
K=(J-2)*(N-2)+I-l 

UP=0.5*(PSI(U+l)-PSI(U-l)+PSI(I+U-i-l)-PSia+U-l)) 

UM=0.5*(PSI(I4+1)-PSI(I4-1)+PSI(MJ+1)-PSI(I-14-1)) 

VP=-0 J*(PSI(I+ 1 4)-PSI(I-l J)+PSI(I+1J+1)-PSI(I-1 4+1)) 

vM=-o.5*(Psi(i+u)-psia-u)+psia+U‘i)-psi(i-u-i)) 

A1(K)=-PE*((UM+ABS(UM))/(2.0*DX**2)+(L0/(PE*DX**2))) 

TEMP=(UP+ABS(UP)-UM+ABS(UM))/(2.0*DX**2) 
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B 1 (K)=PE*((2.0/(PE*DT))+TEMP+(2.0/(PE*DX**2))) 

C 1 (K)=PE*((UP- ABS(UP))/(2.0*DX**2)-( 1 .0/(PE*DX* *2))) 

TEMPI =(VP-ABS(VP))*TH(I J+l )+<VP+ABS(VP)-VM+ABS(VM))*TH(U) 
TEMP 1 =TEMP 1 -(VM+ AB S( VM)) *TH(I,J- 1 ) 

TEMP 1 =(TEMP 1 )/(2.0*DX* *2) 

TEMP3=TH(IJ+ 1 )+TH(I J- 1 )-2.0*TH(I J) 

TEMP3=( 1 .0/(PE*DX* *2))*TEMP3 
D 1 (K)=PE*((2.0*TH(1 J)/(PE*DT))-TEMP 1 +TEMP3) 

IF G.EQ.2) THEN 
A1(K)=0.0 
ENDIF 

IF (LEQ.(N-1)) THEN 
D 1 (K)=D 1 (K)-C 1 (K)*TH(I+ 1 J) 

C1(K)=0.0 

ENDIF 

85 CONTINUE r 

80 CONTINUE 

C 

C SOLVE FOR TEMPERATURES AT HALF TIME STEP 

C 

CALL TRID(NUM f AlJBLCl,DlJEl) 

DO 8U=2,M-1 
DO 86,I=2,N-1 
K=(J-2)*(N-2)+I-l 
TH(IJ)=E1(K) 

86 CONTINUE 

81 CONTINUE 

C 

C NOW MAKE THE Y -DIRECTION IMPLICIT AND SOLVE FOR TEMPERATURE 

C 

DO 82,I=2,N-1 
DO 87J=2 f M-l 
K=(I-2)*(M-2)+J - 1 

UP=0.5*(PSI(U+1)-PSI(U-1)+PSI(I+1J+1)-PSI(I+U-1» 
UM=0.5*(PSI(U+1)-PSI(U-1)+PSI(I-U+1)-PSI(I-U-1)) 
VP=-0.5*(PSI(I+U)-PSI(I-U)-t-PSI(I+l T J+l)-PSI(MJ+l)) 
VM=-0.5*(PSI(I+1 J)-PSI(I- 1 J)+PSI(I+1 J-l )-PSI(I-U- 1 )) 
Al(K)=-PE*((VM+ABS(VM))/(2.0*DX**2)+(1.0/(PE*DX**2))) 
TEMP=(VP+ABS(VP)-VM+ABS(VM»/(2.0*DX**2) 

B 1 (K)=PE*((2.0/(PE*DT))+TEMP+(2.0/(PE*DX**2))) 
C1(K)=PE*((VP-ABS(VP))/(2.0*DX**2)-(L0/(PE*DX**2))) 

TEMPI =(UP- ABS(UP))*TH(I+ 1 , 1 )+(UP+ABS(UP)-UM+AB S(UM))*TH(IJ) 
TEMP1=TEMP1-(UM+ABS(UM))*TH(I-1 J) 
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TEMPI =(TEMP1 )/(2.0*DX**2) 

TEMP3=TH(I+ 1 J)+TH(I- 1 r J)-2.0*TH(I r J) 

TEMP3=( 1 .0/(PE*DX**2))*TEMP3 
D1(K)=PE*((2.0*TH(U)/(PE*DT))-TEMP1+TEMP3) 

IF (J.EQ.2) THEN 

IF ((I.GE.N 1 ). AND.(I.LE.N2)) THEN 
B 1 (K)=B 1 (K)+(4.0/3 .0)* A 1 (K) 

C1(K)=C1(K)-(1.0/3.0)*A1(K) 

D 1 (K)=D 1 (K)- A 1 (K)*(2.0*DX/3 .0) 

A1(K)=0.0 

ELSE 

Bl(K)=Bl(K)+(4.0/3.0)*Al(K) 

C 1 (K)=C 1 (K)-( 1 .0/3 .0)* A 1 (K) 

A1(K)=0.0 
ENDIF 
END IF 

IF (J.EQ.(M-1)) THEN 

B 1 (K)=B 1 (K)+(4.0/3 .0)*C 1 (K) 

A 1 (K)= A 1 (K)-( 1 .0/3 .0)*C 1 (K) 

C1(K)=0.0 

ENDIF 

87 CONTINUE 
82 CONTINUE 

CALL TRID(NUM,A1,B1,C1,D1,E1) 

DO 88,I=2,N-1 
DO 89J=2,M-1 
K=(I-2)*(M-2)+J - 1 
TH(IJ)=E1(K) 

89 CONTINUE 

88 CONTINUE 
C 

C USE ALTERNATING IMPLICIT DIRECTION SCHEME TO SOLVE FOR 
VORTIOTY 

C AND START WITH THE X-DIRECTION BEING IMPUOT 
C 

DO 90,J=2,M-1 
DO 1004=2^-1 
K=(J-2)*(N-2)+I-l 

up=0.5*(psi(u+i)-psia4-i)+psia+u+i)-psia+u-i)) 

uM=o^*(Psi(u+i)-psia4-D+psia-u+i)-psia-u-i)) 

vp=-o^*(Psia+u)-psia-u)+psi(i+u+i)-psia-u+i)) 

VM=-0.5*(PSI(I+1 4)-PSia-U)+PSI(I+l 4 -l)-PSI(I-U- 1 » 
A1(K)=-PE*((UM+ABS(UM))/(2.0*DX**2)+(1.0/(RE*DX**2))) 
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non 


TEMP=(UP+ AB S(UP)-UM+AB S(UM))/(2.0*DX* *2) 

B 1 (K)=PE*((2.0/(PE*DT))+TEMP-K2.0/(RE*DX**2))) 
C1(K)=PE*((UP-ABS(UP))/(2.0*DX**2)-(1.0/(RE*DX**2))) 
TEMP1=(VP-ABS(VP))*W(I,J+1)+(VP+ABS(VP)-VM+ABS(VM))*W(I,J) 
TEMPI =TEMP1 -(VM+ABS(VM))*W(U-1 ) 

TEMPI =(TEMP1)/(2.0*DX**2) 

TEMP2=(TH(I+ 1 r J)-TH(I- 1 ,J))*SIN( ALPHA*PI/1 80.0) 
TEMP2=TEMP2-(TH(U+ 1 )-TH(U-l ))*COS(ALPHA*PI/l 80.0) 
TEMP2=-(RI/(2*DX))*TEMP2 
TEMP3=W(U+1 )+W(I,J- 1 )-2.0*W(U) 

TEMP3=(1 .0/(RE*DX**2))*TEMP3 

D1(K)=PE*((2.0*W(I,J)/(PE*DT))-TEMP1+TEMP2+TEMP3) 

IF (LEQ.2) THEN 

D1(K)=D1(K)-A1(K)*W(I-U) 

A1(K)=0.0 
END IF 

IF (I.EQ.(N-1)) THEN 
D1(K)=D1(K)-C1(K)*W(I+U) 

Cl (K)=0.0 
END IF 

100 CONTINUE 

90 CONTINUE 

CALL TRID(NUM,A1,B1,C1,D1,E1) 

DO 72,J=2,M-1 
DO 73,I=2,N-1 
K=(J-2)*(N-2)+I-l 
W(U)=E1(K) 

73 CONTINUE 

72 CONTINUE 

MAKE THE Y -DIRECTION IMPLICIT NOW AND SOLVE FOR VORTICITY 

DO 91,I=2,N-1 
DO 10U=2,M-1 
K=(I-2)*(M-2)+M 

UP=0.5*(PSI(U+1)-PSI(U-1)+PSI(I+U+1)-PSI(I+1J-1» 

uM=o.5*(Psi(u+i)-psiaj-i)+psia-u+i)-psia-u-i)) 

VP=-0.5*(PSI(I+U)-PSI(I-l,J)+PSia+U+l)-PSI(I-U+l» 
VM=-0.5*(PSI(I+U)-PSia-U)+PSI(I+U-l)-PSI(I-lJ-l)) 
A1(K)=-PE*((VM+ABS(VM))/(2.0*DX**2)+(1.0/(RE*DX**2))) 

TEMP=( VP+ AB S(VP)- VM+ AB S(VM))/(2.0*DX**2) 
B1(K)=PE*((2.0/(PE*DT))+TEMP-K2.0/(RE*DX**2))) 
C1(K)=PE*((VP-ABS(VP»/(2.0*DX**2H1.0/(RE*DX**2») 
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TEMP1=(UP-ABS(UP))*W(I+U)-KUP+ABS(UP)-UM+ABS(UM))*W(U) 
TEMP1=TEMP1 -(UM+ABS(UM))*W(I-1 J) 

TEMPI =<TEMP1 )/(2.0*DX**2) 
TEMP2=(TH(I+U)-TH(I-U))*SIN(ALPHA) 
TEMP2=TEMP2-(TH(I,J+l)-TH(IJ-l))*COS( ALPHA) 
TEMP2=-(RI/(2*DX))*TEMP2 
TEMP3=W(I+1J)+W(I-U)-2.0*W(U) 

TEMP3=(1 .0/(RE*DX**2))*TEMP3 
D1(K)=PE*((2.0*W(U)/(PE*DT))-TEMP1+TEMP2+TEMP3) 

IF (J.EQ.2) THEN 

D1(K)=D1(K)-A1(K)*W(U-1) 

A1(K)=0.0 
END IF 

EF (J.EQ.(M-l)) THEN 
D1(K)=D1(K)-C1(K)*W(U+1) 

C1(K)=0.0 
END IF 

101 CONTINUE 

91 CONTINUE 

CALL TRID(NUM,A1,B1,C1,D1,E1) 

DO 92,1=24^-1 
DO 93J=2,M-1 
K=(I-2)*(M-2)+J- 1 
W(IJ)=E1(K) 

93 CONTINUE 

92 CONTINUE 

C 

C SOLVE FOR STREAM FUNCTION THROUGH ITERATION 

C 

95 DDFF2=0.0001 

DIF=0.001 
DO 102,I=2,N-1 
DO 103J=2,M-1 
PSIOLD(U)=PSI(U) 

TEMP=PSI(I+U)+PSI(I-1J)+PSI(U+1 )+PSI(U- 1 ) 

PSI(I J)=0.25 *((W (I r J)*DX**2)+TEMP) 

DEFFl=ABS(PSI(IJ)-PSIOLD(IJ)) 

IF(DIFF1 .GT.DIFF2) THEN 
DEFF2=DIFF1 
DEF=DIFF2/PSI(I r I) 

END IF 

103 CONTINUE 

102 CONTINUE 
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I=N 

DO 104J=2,M-1 
PSIOLD(IJ)=PSI(IJ) 

PSI(I,J)=0.5*(5*PSI(I-U)-4*PSI(I-2J)+PSI(I-3 t J» 

DIFFl=ABS(PSI(U)-PSIOLD(U)) 

IF (DIFF1 .GT.DIFF2) THEN 
DIFF2=DIFF1 
DIF=DIFF2/PSI(I r J) 

ENDEF 

104 CONTINUE 
C 

C ITERATE UNTIL STREAM FUNCTION CHANGE IS LESS THAN 
C HALF A PERCENT 
C 

IF (DIF.GT.0.005) THEN 
GOTO 95 
END IF 
C 

C CHECK IF SYSTEM HAS REACHED STEADY STATE 

C 

IF (MODE.EQ.2) THEN 
DO 304J=1,M 
DO 301,I=2,N 

IF (TH(U).GT.O.OOOOOl) THEN 
DIFF1 =ABS(TH(U)-TOLD(IJ)) 

DIFF1=DIFF1/TH(IJ) 

IF (DIFF1.GT.(10.0*DT)) THEN 
DO 302J1=1,M 
DO 303,11=2^ 

TOLDaUl)=THaUl) 

303 CONTINUE 

302 CONTINUE 

GOTO 25 
ENDIF 
END IF 

301 CONTINUE 

304 CONTINUE 
GOTO 125 

ENDIF 

C 

GOTO 25 
C 

C OUTPUT RESULTS 
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c 

125 WRITER, 260)RI,RE,PE, TIME 

DO 1 10J=M,1,-NSTEP 

WRITE(*,270XPSI(U),I=l t N,NSTEP) . 

110 CONTINUE 

DO 111 J=1,M,NSTEP 

WRITE(5,270)(PSI(I J),I=1 .N.NSTEP) 

111 CONTINUE 
WRITE(*,201) 

DO 1 20J=M, 1 ,-NSTEP 

WRITE(*,280)(TH(I,J),I=1 ,N,NSTEP) 

120 CONTINUE 
J=1 

NUT=0.0 

DO 113,K=N1+1,N2 

NUX=((K-N 1 )*DX)/TH(K J) 

NUT=NUT+NUX 
113 CONTINUE 

NU=NUT/(N2-N1) 

PRINT* .'AVERAGE NU = \NU 
DO 12U=1,M,NSTEP 

WRITE(4,280)(TH(U),I=1,N,NSTEP) 

121 CONTINUE 
PAUSE 

close(4) 

close(5) 

C 

C 

200 FORMATC ENTER CHANNEL LENGTH, HEIGTH, ENTRY LENGTH IN CM ) 

201 FORMATC ’) 

202 FORMATC CHOOSE EITHER 1 OR 2: ') 

203 FORMATC 1. ENTER DIMENSIONLESS PARAMETERS (RI,RE,PE)') 

204 FORMATC 2. ENTER DIMENSIONAL QUANTmES ) 

205 FORMATC ENTER N,M,N1 .N2.NSTEP*) 

208 FORMATC ENTER ANGLE OF INCLINATION WITH VERTICAL ) 

209 FORMATC ENTER FLOW TYPE: 1) SLUG FLOW 2) POUISELLE FLOW) 

210 FORMATC ENTER INLET TEMP(K),BETA (1/K), HEAT FLUX (W/cm A 2)') 

211 FORMATC ENTER DTJFINAL TIME.PRINT TIME (1 S = .0001)') 

212 FORMATC ENTER MODE: 1) TRANSIENT) 

213 FORMATC 2) STEADY-STATE’) 

214 FORMATC ENTER HEATER LENGTH, DX IN CM*) 

215 FORMATC RI= \F9.2,’ RE= ’,F9.1,’ PE= 'JF9. 1) 

220 FORMATC ENTER INLET VELOCITY (cm/s), GRAVITY (g) ’) 
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230 FORMATO ENTER THERMAL CONDUCT (W/(m K»,DIFFUSITY (m A 2/s)') 
240 FORMATC' ENTER KINEMATIC VISCOSITY (m A 2/s)’> 

250 FORMATC ENTER RI,REJ>E ) 

254 FORMATC SOLUTION DOES NOT CONVERG! ! !’) 

256 FORMATC N= ,13,’ Nl= ’,13/ N2= ,13,’ M= \I3,’STEP= ,13) 

257 FORMAT(I3,I3) 

259 FORMAT(l 00F6.2) 

260 FORMATC RIs'JFM/ RE^J^.l,’ PE=’,F7.3,' TIME=’,F7.3) 

261 FORMATC TIME = \F9.4) 

270 FORMAT(100F7.3) 

280 FORMATO 00F6.3) 

290 FORMATC’ CONTINUE? (l=NO)’) 

295 FORMATC ISOX DATA FILE:STREAM LINES’) 

300 FORMATC ISOX DATA FILE: TEMP LINES’) 

END 

C 

C TRIDIAGONAL MATRIX SOLVER 
C 

SUBROUTINE TRID(N,D1,D2,D3,R,F) 

REAL*8 D 1 ( 1 0000),D2( 1 0000), D3( 1 0000),R( 1 0000) 

REAL*8 F(10000) 

INTEGER N,I 
DO 15,1=2^ 

D2(I)=D2(I)-(D3(I- 1 )*D1 (I)/D2(I- 1 )) 

R(I)=R(I)-(D 1 (I)*R(I- 1 )/D2(I-l )) 

15 CONTINUE 
C 

C USE BACKWARD SUBSTITUTION TO SOLVE FOR F(X) 

C 

F(N)=R(N)/D2(N) 

DO 25,I=N-1,1,-1 

F(I)=(R(I)-D3(I)*F(I+1))/D2(I) 

25 CONTINUE 
RETURN 
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THE VELOCITY POTENTIAL OF TWO SPHERES MOVING 
PERPENDICULARLY TO THE LINE JOINING THEIR CENTERS AND 
THE ASSOCIATED LIFT FORCE 

I. Introduction 

Vapor bubble* generated on a flat heating surface in forced convection boiling, 
aside from being subjected to the influences of forces involved in pool boiling conditions, 
are also subjected to the pressure gradient in the bulk flow field and to the lift produced due 
to the pressure distribution caused by the velocity redistribution. 

Drew (1987) used a numerical method to compute the lift force for small spheres 
moving in a shearing fluid near a wall, but the result is not convenient to use analytically 
and does not include the condition when the sphere touches the wall. 

An alternative to the numerical solution is a potential flow solution involving the 
use of a potential flow model to approximate the real situation when a bubble is on or near a 
wall with bulk liquid flow parallel to the wall in forced convection boiling. The simplest 
potential flow model which best approximates the real situation is that of a spherical bubble 
attaches to or near a wall, with an ideal fluid sweeping over this infinite plane. Since this 
potential flow model is equivalent to two spheres in contact or separated, and moving 
perpendicularly to the line joining their centers in an infinite ideal fluid, it is natural that 
the solution should be pursued along this line. A number of classical solutions to this 
problem exist. 

The solution of two spheres moving in a perfect fluid was first attempted by Stokes 
(1843), then by Herr Bjerknes (1863), and Hicks (1879). However, it was not until 1887 
that Basset(1887) and Herman(1887) solved the problem for two spheres moving 
perpendicularly to the line joining their centers. Even though they were the first ones to use 
the first order associated Legendre polynomials to solve the problem(Basset used the 
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integral expression and Herman used the derivative expression of the polynomial), they 
did not solve the problem completely; instead, they used only the first degree of the 
polynomial in the final expression for the velocity potential, neglecting all higher degree 
terms, most likely due to the lack of computational means. This gave the impression that 
the method was good only for approximate solutions, and much effort was subsequently 
devoted to seeking more accurate solutions. The problem was essentially solved by Basset 
(1887). Miloh (1977) used the proposition proposed and proved by Basset(1887), which later 
appeared in Hobson(1931). The solution of Miloh for the subject problem is identical to 
that obtained by using Basset's method directly. The only difference is in the solution 
technique: in the formulation of Miloh the potential coefficients are solved from an 
infinite set of equations, while in the formulation of Basset the potential coefficients are 
solved by adding a set of infinite summation series. It will be demonstrated in the 
following sections that the two formulations are equivalent Basset was interested in the 
kinetic energy associated with the motion of the spheres, so no results were given relative to 
the motion interaction forces between the two spheres. Miloh did give this force, but with 
insufficient accuracy. The solution was based on solving an infinite set of equations, and 
the residual normal velocity produced by the velocity potential on the spherical surface can 
be a check on the accuracy of the solution. This check was not conducted by Miloh; it was 
simply assumed that the solution satisfied the prespecified boundary conditions, which 
was not the case. 

Miloh had the advantage of knowing the form of the potential function, and 
therefore was able to set out to solve the coefficients in that function directly; in Basset's 
time the potential function was completely unknown, and instead of using the geometric 
imaging method of single dipoles, as was done by prior workers, he was the first person to 
successfully use the function imaging method to lead to the form of the potential function 
in the spherical coordinate system. 
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II. The construction of the velocity potential 


The following is an extension of the work of Basset The configuration and 
coordinate systems are shown in Fig. 1. Two spheres of radii r bl and r b2 are moving in 

an ideal stationary liquid with respective uniform velocities Uj and U 2 in the negative X 
direction. Although the generality of Uj * U 2 will be retained for some time, the particular 
case of interest is where Uj * U 2 , r bl « r M , and the X direction is perpendicular to the line 
joining the centers of the spheres as shown. For clarity, the notation r bl and r b2 will be 
retained, although it » understood that r bl = r b2> 



Fig. 1 . Two spheres moving perpendicularly to the line joining their 
centers 

Spherical coordinates r, 0, 4. and r‘ ,6', 4 are fixed at the centers of the spheres r^ 
and r b2 at Oj and 0 2 , respectively. The Cartesian coordinate system X,Y, Z has the 
origin at the mid-point between the two spheres, with the X axis perpendicular to this line 
and the angle 4 measured in the X-Y plane. 

Lengendre polynomials with degree n and order m, whose origins are at Oj and 
0 2 , respectively, are defined by the following relations: 
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Pjj 1 = M(ain 0 P J" ( 4 + Vn 2 -1 cosa) n ' m (sina) 2,n da 

x 

P’“ *M(sine7 n J ( jx ' + Vp' 2 -1 cosa^sina) 2 ® da (1) 


or 


(sina) 2m da 


Pj 1 = M(sin0) m f 

W 0i+ Vn 2 -1 

ff 

P'“ = M(8in0') n > J p=i 

w Oi' + Vn’ 2 -! 


cosa) n+ ® +1 
(sincQ ^da 

CO BCx) n+m '*'l 


.(!*) 


where 


M = 


(n+m)! 


(n-m)!l«3»...(2m-l)x 


4 = cos 6 

|X* = COB 6' 

Basset was the first one who proposed and proved the following propositions, which relate 
the associated Legendre polynomials P^ and P*^: 

p n r® .(n+m)! m (n+m+1)! jr m 

r *n+l "(n.m)!c n+m+1 2m! * (2m+l)! c p m+l + "* 


k (n+m+k)! 
l ” (2m+k)! 


S* Ck ♦•••> 


(2) 


when r<c, and 

(_l)n-m p® 


jji+l 


r' m (ntm)! 4 n (n+m+1)! r' ,m 

(n-m)! c n+m +l 2m! p nt + ( 2m+l)! c p m+l + ~ 


(n+ro+k)! jm 
+ (2m+k)! 'T^ P ®+k 


+...) 


(3) 


when r’cc 

Since velocity potentials are linear, the total velocity potential 4> at any points 
outside the spheres is the sum of the potentials <&jCosq> produced by the sphere of radius 
centered at Oj and d^coscp produced by the sphere of radius r^ centered at O2. For the sake 
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of simplicity, here and $2 are used to denote the amplitude of the two velocity potentials, 
i.e. 


(C*! + $ 2 ) COS 9 (4) 

Once «t>i is known, can he determined directly from geometric symmetry. <X>j 
must be constructed in such a manner that it satisfies the boundary conditions on the 
surfaces of the two spheres: 


1 

d*i 

-T-i*. = 0 

or r " T b2 


(5) 

( 6 ) 


Beginning with the well-known potential for a single sphere moving in an infinite 


quiescent pool of incompressible ideal liquid: 
3 

* TT rbl 1 


.(7) 


which satisfies the boundary condition given by Eq. (5). Equation (3) is used to transform 
Eq. (7) into the velocity potential in the vicinity of the sphere centered at O 2 : 

<«> 

2c J j«i c l 1 

Eq. (8), however, does not satisfy the boundary condition Eq. (6). To satisfy Eq. (6), Eq. (9) 

is added to Eq. (8): 

3 s» 2i+! t i 

U| c 1 -' (> 

G >1 in the vicinity of the sphere O 2 is now given by: 

3 m 2i+l 

♦l»Ui^2(r* + (10) 

2c 3 " ,+1 r' ,+1 c* 1 
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Although Eq. (10) now will satisfy the boundary condition Eq. (6), it will no longer 
satisfy boundary condition Eq. (5). To demonstrate this, using Eq. (2), Eq. (9) may be 
transformed into: 

'?,Ui 


3 TT - 2i+l - 

. V!iy f j ^2 1 V (-pHq+j)! ri 1 

’ 2c 3 M i+1 (i-l)!c i+2 ^ <*+!)! ci- 1 P i 

2 m *2 2 itl m 

f » U| V , M V (-pl'^i-fj)! r> l, 

* <i+l»!c Ii+1 P 0+»l c)' 1 ’ 

r bl Ul y y (.pj- 1 ri 1 i 2 (i+j)! r b2 
= 2c 3 ffg 1 (j+1)! cJ' 1 ^ (i+D! c 2i+1 


.( 11 ) 


upon adding Eq. (11) to Eq. (7), the velocity potential <X>i in the vicinity of sphere Oj is now 
3 3 .. — — 2i+l 

y 1 T b> Ul V V (-1V' 1 rJ l i 2 (i+j)! r b2 
° 1_Ul 2r 2?1+ 2c 3 jTi' -T? (i+Dicj* 1 ** (i+1)! c 2i+1 (12) 


Obviously, the sum of the two terms on the right hand side of Eq. (12) cannot satisfy 


Eq. (5). 


To satisfy the boundary condition given by Eq. (5), Eq. (13) is added to Eq. (12), 


3 TT - - 2j+l 1 2i+l 

r bi u i y y (.pi-ij r bi p; i2 (i+j)! n>2 

2c 3 jTi'lrf (j+l>(j + D! ei' 1 ri +1 G+1M c 2i+1 


(13) 


The velocity potential in the vicinity of the sphere at Oj is now 
1 _3 „ - ~ 2j+l 1 


>l*Vl 


bl 


1 

P, + 


T bi u i y V 


i -1 


2i+l 


2T 2 1 2c 3 


. , , M- 
U £ 0+H! 


(ri + r 


J 3 i^i+j)! r b2 


j+lji+lV 1 (i+1)! c 2i+1 


(14) 

By summing Eqs. (7), (9) and (13), the potential function produced by the motion of 
sphere Oj, valid at all locations and still satisfying Eq. (5) but not Eq. (6), is: 


143 



3 3 ~ 2i+l ,1 

* „ f bl 1 TT r bl V i f b2 P i 

®i = U, — “ p, + U, A — r + 

1 1 2? 1 1 2c 3 1+1 c M r' i+1 


.(15) 


3 T t •• •• 2j+l 1 2i+l 

M U| V y (.l)Hj r bl P) j2(j+j)! r b2 

+ 2c3 (j+l)(j+D! cH jj+1 (i+1)! c 2i+1 

The above procedure can now be expanded, alternately satisfying Eqs. (5) and (6), 

resulting in Eq. (16), the general potential function for sphere Oj valid at all locations. 

The convergence of this function is shown by the convergence of the normal velocity, i.e., 

its derivative, at the sphere surface in section IV below. 

3 3 «• 2i+l ,i 

- it I^l_ 1 ,y i f b2 P i 

1 = U 1 2r 2 P l + U 1 2c 3 { ^-i + l c i-l r .i + l + 

- - 2j+l 1 2i+l 

*yy — i jbi_^faij)!jb2_ 

W M (j+lXj+U! ci-l ,1*1 (i+1)! ,21+1 * 

- » " 2k+l ,1 2j+l 2i+l 

yyy k ^ i^o+-v>! *i ga+j); r i>2 

*&i & M (k+l)(k+l)! ck-l ,.k+l ( (j+l)!)2 c 2j+l (i+1)! c 2i+l + 

~ ~ 21+1 1 2k+l 2j+l 

+y yyy <1+1 1 bi * po+m rti , 

w Ef j^ITT 0+lK}+D! e l-l r'*‘ ((k+l)!) 2 c21+l «j+l)!)V’*> 

2i+l 

i 2 (i-*-j)! *2 1 

(i+1)! c 2i+l + ** (16) 

By symmetry <b 2 can be written directly: 

3 3 21+1 1 

^ ,, r b2 ,1 „ r b2 ,V + .+!.i i *bl *1 

^2*^2 , P,+ U 2 { (-1) 1 1 t-t TT ♦ 

2r' 2 1 2c 3 1+1 c* 1 r* +1 

- - 2j+l ,1 2i+l 

y y j b2 Pj j2(j»j); r bl 

+ jw Irf <) +1 ><j + l> ! ei* 1 r*j +1 <i+D! c 2i+1 * 

- - - 2k+l 1 2j+l 2i+l 

yyy ,.i*> * m Tb2 rbi , 

&1 ft M (k+l)(k+l)! ^-1 ,k+l ,(,,1)1)2 c 2j+l (i+1)! c 2i+l 

— 21+1 ,1 2k+l 2j+l 

y y y y i ^2 >i h^h+Di'bi j 2 a+h)! r b2 

£f jff M (1+1>(1+1)! C 1 * 1 r' ,+1 ((k+l)!) 2 c 2k+1 ((j+1)!) 2 c 2 -i +1 * 
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2 i+l 

i 2 (i+j)! I * )1 


(i+1)! c 2i+l + "' 


Substituting Eqs. (16) and (17) into (4), and collecting terms that have the same p* 


^ Uj+'costp 


where 4' is the amplitude for the potential of unit velocity in the direction of Uj, given by 


An's and B n 's in (19) are given by the following infinite summation series 
3 n-1 3 •* n-1 2i+l 

^2 b2 n bl bl V n r bl i 2 (i+n)! r b2 
A " = U, 2c 3 n+1 c n-l + n+1 c n-l (n+l)!(i+l)! c 2i+l + 

3 •• •* n-1 2j+l 2i+l 

U2 b2 y y n r bl j 2 (j+n)! r b2 i 2 (i+j)! r bl 

+ U i 2c 3 n+1 c n Un+l)!(j+l)! c 2j+t(j+l)!(i+l)! c ^ 

_3 ~ — n -l 2k+l 

— T y Y n Tbl k 2 (k+n)! r b2 j 2 (j+k)! 

+ 2c r* T* n+1 c n-l (n+l)!(k+l)! c 2k+l (k+l)!(j+l)! 


k*l j*l i*l 
2i+l 

i 2 (i+j)! r b2 
'(j+l)!(i+l)! c 2i+l 
3 n-1 


n-1 2k+l 2j+l 

n T bl k 2 (k+n)! r b2 j 2 (j+k)! T b* 
n+1 c n-l (n+l)!(k+l)« c 2k+l (k+l)!(j+l)! c 2j+l 


- 3 •• n-1 2i+l 

„ bl n b2 U2 r b2 y n r b2 j 2 (i+n)! r bl 

n ~ 2 c* n+1 c”* 1 + U| n+1 c”* 1 (n+l)!(i+l)! c 2i+l + 

3 — ~ n-1 2j+l 2i+l 

^iV y n r b2 j 2 (j+n)? r bl i 2 (i+j)! r b2 
2c 3 ^ n+1 c"" 1 (n+l)!(j+l)! c 2j+l (j+l)!(i+l)! c 2i+l + 

3 •• — •• n-1 2k+l 2j+l 

U 2 b2 y y y n r b2 k 2 (k+n)! r bl j 2 G+k)! r b2 
+ U l2c 3 n+ l c® 1 (n+l)!(k+l)! c 2k+l (k+l)!(j+l)! c 2j+l 


3 j» •• n-1 2j+l 2i+l 

^biy y n r b2 j 2 (j+n)? r bl i 2 (i+j)? r b2 

2c 3 ^ n+1 c"" 1 (n+l)!(j+l)! c 2j+l (j+l)!(i+l)! c 2i+l + 

3 m mm m n-1 2k+l 


Ul j-1 t-1 
2i+l 

i 2 (i-+j)! r bl 
*g+D!(i+i)! c 2i+i + - 


n-1 3 3 •• 2 

n r bl ^^b2 r biy j 2 (j+ n )! r b 

^“n+l c n-l *U^ 2c 3 + 2 <?^f (n+l)!(i+l)! c 2 
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3 — 2j+l 2i+l 

J^2__^y j 2 (j+n)! r b2 i 2 (i+j)l T *>1 

+ U 1 2c 3 ff (n+l)!(j+l)! c 2j+l (j+l)!(i+l)! c 2i+l + 

3 2k+l 2j+l 

— y y y ^*±£>1 14,2 

+ 2c 3 W H M (n+1)!(k+1)! C^ 1 (k+DKj+D! c 2 i +1 * 
2i+l 

, i 2 (i+j)! r b2 

(j+l)!(i+l)! c 2i+l + 

n-1 3 3- 2i+l 

n b2 *bl U 2 r b2 V i 2 (i+n)l r bl 
““'n+l c n-l 2 C 3 + U, 2 c 3 ^ (n+l)!(i+l)! e 2i + l + 

3 - - 2j+l 21+1 

_Liy y j 2 (i+n)l r bl i 2 (i+j)l r *>2 

+ 2c 3 ^ -ff (n+l)!(j+l)! c 2j+l (j+l)!(i+l)! c 2i+l + 

3 2k+l 2j+l 

U 2 b2 V V y k 2 (k+n)! r bl j 2 (j+k)! r *>2 

+ U l2c 3 uiH ^ (n+1)Kk+1)! c 2k+1 (k+l)!(i+l)! c 2j+l 

2i+l 

, ** , 

(j+l)!(i+l)! c2i+i 


( 20 ’) 


( 21 ') 


In (20) and (21), Uj and U 2 are the parallel velocity components of sphere 1 and 
sphere 2, respectively. They may be identically the same or different in either or both 
direction and magnitude, and the positive direction of Uj is in the negative X direction. In 
the case that either the direction or magnitude or both of one velocity is different from that of 
the other, the above result is correct only at the moment when the two spheres are moving in 
such a way that these two velocity components are perpendicular to the line joining their 
centers. 


Ajj's and B n 's in (20) and (21) can be evaluated one at a time. However from (20) and 
(21) it can be seen that the result of a lower summation order term in Ag could be used in the 
evaluation of a next higher summation order term in B„ and vice versa. It is obvious that 
it would be more advantageous to evaluate these simultaneously in groups (with the upper 
limit of N prespecified) in order to reduce computational effort. These have been computed 
below for two identical spheres, with U 2 , and for the cases of various center distances. 
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For the case of two identical spheres touching and moving with the same velocity, 
U 1 =U 2 , r bl = r k 


b2 


r and c=2r Eq. (19) reduces to: 
b b* 


r 3 r 3 
.. b 1 r b 

V 


2r* 


2r-‘ 


•;*X a , c 2 ( 


n«l 


1 ,1 

Pa 

jj»+l + r 'B+l 


.( 22 ) 


where 


. 1 . 1 n 1 n y i 2 (i+n)! 1 

2* ' 2 n ' ln+1 + 2 n * ln+1 (n+DKi+DJg 2 ** 1 + 

1 n V V j 2 (j+n)! 1 i 2 (i+j)! 1 

+ 2 n - 1 n+l^j' “ (n+l)!(j+l)! 2 2j+l(j+l)Ki+l)!2 2i+1 + 

as es as 

J__n_V y y k 2 (k+n)! 1 j 2 (j+k)! 1 

+ 2 n * ln+1 j~J < ^ ^ (n+l)!(k+l)! 2 2k+1 (k+l)!(j+l)! 2 2 J +1 * 

. ihi±m , 

(j+l)!(i+l)! 2 2 >+l •**' 


(23) 


or 
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2 n+3 n+1 
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♦X 
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i 2 (i+n)! 1 

(n+ !)!«+!)! 2 2i+1 


y y j 2 (i+n)! 1 i 2 (i+j)! 1 

+ (n+l)!(j+l)! 2 2j+i (j+l)!(i+l)! 2 2i+1 + 

y y y k 2 (k+n)i 1 j 2 (j+k)i 1 i 2 (i+j)! 1 , 

+ " H (n+l)!(k+l)! 2 2k+l (k+l)!(j+l)! 2 2 i +1 (j+DKi+l)! 2 2i+1 


(23’) 


The A„'s in Eq. (22) were evaluated by using Eq. (23). Table 1 gives the first 16 
terms for the velocity potential with different numbers of terms. The Aq's were computed 
for N up to 2000. The relative error for each term in Eq. (23) was set as 10* 17 , and the 
relative error for A^s was set as 10* 1 1 . Eq. (23) was also used to evaluate the A^s for N up 
to 30, with relative error set as 10* 4 . Since each term in Eq. (23) uses the results of the 
previous term no relative error was set for each term in Eq. (23) in this case; instead each 
term is computed for the index from 1 to 30. The results are identical to those given by 
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Miloh (1977) . This reveals why the values given by Miloh were smaller than what is 
obtained as the total number of terms is increased. Even though N was chosen to yield a 
m a xim u m relative error of 110* 4 between successive approximations in the computations 
by Miloh, the total relative error is larger than ±10~*, as can be seen in the values of the 16th 
term for n*30 and n«1000 in Table 1. 

For the case of two identical spheres moving in the same direction with the same 
velocity and not touching, the coefficients are computed by using Eq. (20') for center 
distances C/R^ = 2.000001, 2.000002, 2.000003, 2.000004, 2.000005, 2.000007, 2.00001, 2.000015, 

2.00002, 2.00003, 2.00005, 2.00007, 2.0001, 2.0002, 2.0003, 2.0005, 2.0007, 2.001, 2.0015, 

2.002, 2.003, 2.005, 2.007, 2.01, 2.015, 2.02, 2.02244, 2.03, 2.05, 2.07, 2.10, 2.20, 2.3, 2.4, 2.5, 
2.6, 2.683, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.415, 3.5, 3.683, 4.0, 4.45, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 
8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 1L5, 12.0, 20.0 and 22.0. When the center distance 
between the two spheres is large(C/R^ 22.02), the computation of the coefficients using Eq. 
(20') converges very rapidly. The first 16 terms of the coefficient in the velocity potential 
for a few selected center distances are given in Table 2. The results for the remaining 
eases listed are available, but not given here. 



Table 1. 

The first 16 coefficients 

II 

jf 

5 

<* 

2.000) 


n\N 

30 

100 

200 

500 

1000 

2000 

1 

.04032383 

.04034059 

.04034122 

.04034132 

.04034133 

.04034133 

2 

.02998921 

.03001425 

.03001519 

.03001535 

.03001536 

J03001536 

3 

.01952214 

.01955302 

.01955418 

.01955437 

.01955438 

.01955438 

4 

.01258031 

.01261584 

.01261716 

.01261738 

.01261740 

.01261740 

5 

.00630055 

.00634003 

.00834151 

.00834175 

.00834177 

.00834177 

6 

.00567546 

.00571841 

.00572002 

.00572028 

.00572030 

.00572030 

7 

.00403218 

.00407826 

.00407998 

J00408027 

.00408029 

.00408029 

8 

.00296983 

.00301876 

.00302059 

.00302089 

.00302092 

.00302092 

9 

.00225711 

.00230868 

.00231061 

J00231093 

.00231096 

.00231096 

10 

.00176079 

.00181483 

.00181685 

.00181718 

.00181720 

.00181721 

11 

.00140299 

.00145934 

.00146145 

J00146180 

.00146182 

.00146183 
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12 

.00113708 

.00119562 

.00119781 

.00119817 

.00119819 

.00119820 

13 

.00093425 

.00099487 

.00099714 

.00099751 

.00099754 

.00099754 

14 

.00077610 

.00083872 

.00084106 

.00084145 

.00064148 

.00084148 

15 

.00065050 

.00071504 

.00071745 

.00071785 

.00071788 

.00071789 

16 

.00054916 

.00061556 

.00061804 

.00061845 

.00061848 

.00061848 


Table 2 . The first 16 coefficients A vs dr. 

n b 


nV 

— =2400000 
r b 

1000001 

1000010 

1000100 

1001000 

1010000 

1100000 

1200000 

1 

.04034133 

.04341190 

.04339905 

.04032728 

.04020703 

.03915365 

.03195988 

.02664009 

2 

.03001536 

.03001519 

.03001369 

.02999903 

.02986126 

.02869903 

.02153822 

.01679388 

3 

.01955438 

.01955422 

.01955274 

.01953837 

.01940547 

.01833209 

.01245318 

.00902192 

4 

.01261740 

.01261725 

.01261584 

.01260231 

.01247927 

.01153133 

.00696264 

.00464970 

5 

.00834177 

.00834162 

.00834030 

.00832765 

.00821460 

.00738535 

.00388762 

.00237215 

6 

.00572030 

.00572016 

.00571892 

.00570706 

.00560291 

.00487594 

.00219639 

.00121337 

7 

.00408029 

.00408016 

.00407898 

.00406780 

.00397134 

.00333046 

.00126340 

.00062634 

8 

.00302092 

.00302080 

.00301967 

.00300908 

.00291927 

.00235087 

.00074188 

.00032747 

9 

.00231096 

.00231084 

.00230976 

.00229968 

.00221566 

.00170883 

.00044490 

.00017376 

10 

.00181721 

.00181709 

.00181606 

.00180643 

.00172752 

.00127349 

.00027221 

.00009364 

11 

.00146183 

.00146171 

.00146071 

.00145150 

.00137714 

.00096886 

.00016959 

.00005125 

12 

.00119820 

.00119809 

.00119712 

.00118827 

.00111801 

.00074966 

.00010734 

.00002846 

13 

.00099754 

.00099744 

.00099650 

.00098798 

.00092144 

.00058816 

.00006886 

.00001602 

14 

.00084148 

.00084138 

.00084047 

.00083225 

.00076911 

.00046681 

.00004468 

.00000913 

15 

.00071789 

.00071779 

.00071690 

.00070896 

.00064893 

.00037413 

.00002926 

.00000525 

16 

.00061848 

. 00061838 

.00061753 

.00060984 

.00055267 

.00030236 

.00001931 

.00000305 


III. The transformation of Equations (20) and (21) into an infinite set of equations 

As was pointed out in the introduction, the formulation of the potential coefficients 
by an infinite series for each coefficient as given by Eqs. (20) and (21) is the same as the 
formulation of Miloh, in which the coefficients are given as an infinite set of equations. 
This will be demonstrated now. 

First rearrange Eq. (19) as: 
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^ n 1 p 1 

a ' * / {(.ly 11 A' r n+2 — — 4-B’ T Jl+2 — ) 

f £-d u ir A ® r bl jj+1 + ®n r b2 . n +l ' 
n«l * * 


where 


An = 2 5 ( n *l) + An 
U 2 

B ’n = 5iT8(n-l) + Bn 


and 6(n-l) it the Kronecker delta function which is one(l) for n*l, and zero(0) otherwise. 
Next rearrange Eq. (25) and substitute X A , 3Lg=-^ , and e = ^ : 

a- I su n ^ A ,t. 3, ... j* V* m^m+n)! 2m+l 

An'I^ 1) (n+l)(n+l)! l2*B (n+1)+ 2 (m+1)! + 

m*l 

as es 

IjjV m 2 (m+n)! .2m+l i 2 (i+m)! .2i+l 

2 ^ (m+1)! ** (m+l)!(i+l)! » + 

m*i i*i 

j, 3 •• « •• 

y yr m 2 (m+n)! .2m+l j 2 (j+m)! .2i+l i 2 (i-»-j)! -2i+l » 

2 , ri (m+D! (m+l)!(j+l)! ** (j+l)!(i+l)! ** + "' 

m*l j*l i*l 

% n- 1 — * m-1 

1 fff „ n , A f e.3 , ... \* (m+n)! .m+2 f m 8 

*2 8(n-l)+ (n+1)(n+1)! (j (n+D- +Zrf ( m .i)! *(m+l)(m+l)! 

111*1 

,3 

i \i i 2 (i+m)! . 2i+l 

[ j (m+1). + (i+1)! ^ + 

i«l 

itf i 2 0* m >' iM !?»> * -m , M1 

2 jMi ml <i +1)! A 0+l)Ki+l)! ^ " J " tZ6) 


Since 


(m+n)l , m+2 


a 3 » «\i t m* - \ +nji . n 

2 Xj (n+l). * j ®(m*l) *b 


substitute Eq. (27) into Eq. (26): 


nX ~ ^ 

A' _ 1 n , A re .,,m+2 (m+n)! V (m+n)! m+2 

A n - 2 8(n * 1) + ( n +l)(n+l)! Ig * 1)X B ( m -l)! ''' (m-l)! 
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, mX B r ^ , ... t,3V i 2 (m+i)! ^ 1 . 

* { (m + lKmTIj! C T (m+1)! + 2 Zrf ( i+ i)! V 

*A y y j 2 (m+j)! . 2j+l i 2 (j ■»•!.)!._ j>*-1 ijj 

+ 2 4i Af 0+1)! X A (j + l)!(i+l)! * * ' in 

J-l 1«1 

.n-1 » »m-l »3 

1 .. . ° A V (m + n)! i™ +2 |?w m n + - r— (m+l)'+ 

“^"'^(n+lXn+l)! ^ (m-1)! *B '2 Km ' l) + (in+lXiB+l)! l 2 ( 

m*l 

e .3 Y l£ls°±i21 ,21+1 V* j 2 (*n+i)! ,2i+l i 2 (j+i)! j2i+l -h 

+ 2*B" (i+l)! X A + 2 A* Ai (j+1)! ** (j+l)Ki+D! ^ ”' J 

i-l J*1 1*1 

(28) 

mXj 1 £ t.iY i 2 (m+i)! .21+1 

B 'm= j 6(m - 1)+ (m+l)(m+l)! { 2 (in+1) - + 2 X *»2rf (i+l)! X A 

VaY V j 2 (m+j)! .2j+l i 2 (j+i)! j2i+l j (29) 

+ 2 (j+1)! X A (j+l)!(i+l)! ® + 

Substitute Eq. (29) into Eq. (28), yields: 

.n-1 — 

A , 1 c,_ ,, " A V n . (*°+ n > ! ^m+2 f3m 

A n = j 5(n ' 1) + (n+l)(n+l)! " Bm (m-1)! 

m*l 

In a similar manner it could be demonstrated that 
nX B (m+n)! m+2 

»’»- I ^^ C.lXn TI)! *'« *A (31) 

m*i 

These are equivalent to the results given by Miloh. 

IV. The evaluation of the residual normal velocity produced by Eq. (19) at the 
surface of the spheres 

The adequacy of the velocity potential constructed above, or others constructed by 
previous researchers, can only be checked by the evaluation of the residual normal 
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velocity, which should be zero on the spherical surface. Since the procedure for evaluating 
the normal velocity on the two spheres are the same, it would be sufficient to evaluate the 
residual normal velocity on either one. The residual normal velocity will be evaluated 
here on the sphere with center at 0}. 


For any given 9 the normal velocity varies with the angle 9, and is a maximum 
when 9=0° or 9*180°, that is, in the plane parallel to the direction of movement of the spheres 
containing both centers. It is sufficient to compute the residual normal velocity for 9=0°, 
with 9 varying from 0° to 180° since the ones for 9=180° are the same in magnitude as those 


for 9=0° and opposite in sign. 

d*' 


v' I =f=-l, 
w r*r bl dr 1 


bl 


.(32) 


or 


nr m, 


bl 


1 u 23 
p i + u ~Tt ( 


r b2 P l 


2r' 2 


•)i„ +X fa+u {(-i^aV* 

bl n-l 


n+2 ,| 

The recursive scheme given by Hobson(1931) ( p290, Eq. 164) 
(n-m) p™(u) = (2n-l)pp” 1 (p) -(n-l+m)p“ 2 (p) 


(33) 


(34) 


together with the initial values 
pj = sin0 

p i*° 

is used to evaluate in Eq. (33). 

Using the definition for p*^ given by Eq. (1) taking m=l and dividing it by r' n+1 , we have: 


.(35) 


.1 


p “ n(n+l) 


r ,n+ l"*r ,2n+1 


% 

sinO’r'^J (p’r'+ V p' 2 -1 r'cosa) n l sin 2 ada (36) 


Substituting the relations in Eq. (37) into Eq. (36), Eq. (38) is obtained. 
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rsinO & r'sinO' 


rp = c + rp 


r' 2 sr I + c 2 + 2cT cos0 
n(n+l)rsin8 


.( 37 ) 


P’i 


.•n+ 1 


itfi^+c 2 +2 c r p) 2 


— j- J (pr+ c + y/n 2 -1 r cosa J”' 1 sin 2 ada 
n+x <f 


.(38) 


The left tide of Eq. (38) is a function of r and 6, and can be readily differentiated 


with respect to r. Performing this operation and evaluating at rer^ yields: 


r n+2 n’ 1 

!(-“ 2-)| . 

3rn+l r . n +l 


n+2 . a 
r b2 sin0 


4-X)n.^.l 

*ii bl 


tonxi* n> 


I|, ' r bi I " l) 


r bl 


r bl 


r bl 


bl 


.(39) 


where 


n(n+l) 


a 

J ( 4 + 7 - +Vp*"- 

<f r bl 


1 cosa) 0 ' 1 sin 2 a da 


.(40) 


It can be proved that I„ can be computed by the following recursive relation 

2 

I »*xr ,<i!n • 1,( ^*> ,) Vl•"« lt x +2 ^7' i)I ‘> , (41) 


'n ‘ bi 


with the initial values 
io=o 

I|*l (42) 

2 n+ — 

Since the evaluation of (1 + + 2“~ \i ) 2 in Eq. (39) is difficult or impossible when n 

*ii Fbl 

is large, Eq.(39) is rearranged as: 
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(43) 


n+2 ,1 

Li.), . 

dr n+1 r . n +l « bl ‘ 


n+2 . 
r b2 Bin0 

: n+2 

r bl 


c 

(-V-l)n- — n-1 


bl 


r bl 


n+1 


K .-rr K ni 1 


'bl 


where 


.(44) 


U + 4 " + 2 ■— 4 )” + 2 

r bl Fbl 


Kq can be computed by the following recursive procedure 

«n T ( 

- bl 


(n+1) (1 + -z- + 2 — n ) 

■ K 1 


'bl bl 


.(45) 


with the initial values of K being given by 
1 


K 


1 


(1 + 4 - + 2^-4 )2 


'bl 


bl 


Ko = 0 (46) 

The residual normal velocity at the sphere surface produced by the second term and 
the series in Eq. (19) are plotted for two identical touching spheres for the cases of 30, 100, 
200, 500, 1000 and 2000 terms in the velocity potential in Figures 2-7 respectively. The 
normal residual velocity computed by using the first 30, 100, 200, 500 and 1000 terms of the 
velocity potential having 2000 terms are plotted in Figs 2b through 6b. 

From Figs 2 to 7 and 2b to 6b it can be seen clearly that as the number of terms in the 
velocity potential increases the amplitude of the residual normal velocity at the sphere 
surface diminishes and the region having higher amplitude narrows. In Fig. 7 the largest 
amplitude is about 2% and the region of high amplitude is about 2 degrees. 

Comparing Figures 2 through 6 with 2b through 6b, one can see that the normal 
residual velocities in the latter group are considerably larger than those in the former 
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group. Therefore, if a velocity potential with fewer terms is desired, the coefficients in that 
potential should be computed specifically for that potential and it should not be constructed 
by simply truncating a potential that has more terms. 

The residual normal velocity at the sphere surface produced by the second term and 
the series in Eq. (19) are plotted for two identical spheres very near to each other, using 2000 
terms in the velocity potential, for various separation distances, in Figures 8 - 13. When 
the center distance between the two spheres is greater than 2.00002, the constructed velocity 
potentials produce a practically zero normal residual velocity, and therefore are not 
plotted. 

It should also be noted from Figs 2 through 13 that the normal velocity as computed 
from Eq. (19) converges and therefore Eq. (19) converges. 


V. The accuracy check of the recursive schemes of Equations (34) and (41) 


It can be seen that in each of the two recursive procedures Equations (34) and (41) and in 
Equation (39) two terms are being subtracted one from the other. If in any step the two terms 
have the same magnitudes and with the first few digits being the same that many 
significant digits would have been lost in that step. Since nearly 2000 such steps were 
performed, it is therefore prudent to ask how accurate the results given by the recursive 
schemes are. 

To verify the value of the associated Lengendre polynomials given by the recursive 

scheme Eq. (34), the following direct scheme was tried: 

a 




n(n+l) 


•in0 J (p. + Vn 2 


1 cosa) n ' 1 sin 2 a da 


n(n-fl) ( n "- 1 \ H a - 3 (p 2 -l)|j+ (n ' 1) 4 V (n ' 4) p n - 5 (n 2 -l) 2 |f - 

<n " 1 (2k)j’ 2k> * i ”" 2K ~ 1 ( ^* 1)K H-^r 2(di) + - ),in0 (2k * n * 1} 

(47) 
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When the order of n is large, some terms in Eq. (47) will have a large magnitude 
and cause a big roundoff error. For example, when n=1000, the largest terms in Eq. (47) 
have the m a gn itude of 10 89 : with double precision only 16 significant digits are available, 
and the roundoff error is of the order of 10^ 3 . Eq. (47) therefore can not be used for large 
orders of n. 

To facilitate numerical computation, the complex integrals of the first order 

Associated Lengendre polynomials as defined by Eq. (3) with m*l are transformed into 

real integrals as given by the following 

x 

pil = n(n±!) sine J (p + "v/p 2 -l coscO^sin 2 *! da 

f — 

- sin 9 j (p 2 sin 2 a + cos 2 a ) 2 (coses + i since )sin 2 a da 


s n(n+l) sin0 


| J «A. 


sin 2 a + cos 2 a ) 2 coses sin 2 a da 


where 


(n-1) Arctan( ^ cosa) = (n-1) Arctan( tand cosa) 


for cosO >0, and 


css(n-lXx +Arctan(^g cosa)) =(n-l)(x+Arctan(tan0 cosa)) 


.(49') 


for cos9<0. 


Special care must be exercised for 6*90° . Examination of Eq. (47) reveals that at 


0*90°: 


when n is even 


* i — * -j 

Simpson's composite rule was first used to evaluate Eq. (48), and the result was not 
satisfactory for higher orders of p*, Romberg's algorithm was then used, and was much 
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superior to Simpson's composite rule. The results of the numerical integration show that 

the recursive scheme Eq. (34) giveB at least 12 accurate digits. 

To verify the accuracy of the result of Eq. (39) with the recursive scheme of Eq. (41), 

Eq. (39) is first rearranged in the following manner so that numerical overflow can be 

avoided during the computation: 

“♦ 2 «•! _ _«+ 2_. 


d, r b2 P; .. nr M sin9 


± 


)l. 


dr'n+1 .n+l w b r , 

r 01 n+2,. c 2 


r n ( 1 + t + 2 t :' 1)2 
r n bl 


c r w 


.(51) 


where 


J n = 


»(i+ 

r bl 


f 

(f r bl 


2 -1 cosa)®* 1 sin 2 a da 


.(52) 


Since J n is a complex integral, J n is transformed to a real integral to facilitate the 


numerical computation: 

J * (p + — l 2 + sin 2 0 cos 2 a na 

{ } 2 (COSO)' 

( l + ( — ) 2 + 2- £ -p) 2 


+ i sinco' ) sin 2 a da 


r bi 


r bl 


(p + — ) 2 + sin 2 © cos 2 a iM 


J ip — r + sin-w cos t i 

{ — — } 

a btr-) 2 * 2 ^ 11)2 


2 cosco' sin 2 a da (53) 


r bl 


r bl 


where 


«■- (n-l) Arctan( ) 


.(54) 


cos© + — 
r bl 

Again, Romberg's Algorithm was used to compute the integral J n and is used in 
Eq. (39). The result shows that the recursive procedure gives at least 12 accurate digits. 
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VL The interacting force on a sphere 


The interacting force on the spheres could be determined by integrating the local 
pressure on the sphere surface around either one of the two spheres, since the force on the two 
would be equal in magnitude and opposite in direction. However, since the basic problem 
is unsteady and the expression for pressure would involve a time dependent term, it would 
be convenient to add a uniform flow in the X direction in the velocity potential as given by 
Eq. (4) and thereby translate the unsteady problem of two spheres moving in the negative X 

direction into a steady problem of two fixed spheres in a uniform flow field. 

= Uj r sin0 cos<p + d> *..(55) 

Since the lift force exerted on a sphere due to the interaction between the sphere and 
a plane is the same as the interaction force between two spheres when the sphere center-to- 
plane distance equals half the center>to-center distance between the two spheres, henceforth 
the term "lift force" will be used indiscriminately for both cases. 

The velocity components in the r, 0 and 9 directions can be determined by 
differentiating Eq.(55). Since the velocity potential given by Eq. (65) is not exact the r 
direction velocity on the sphere may not be zero. However, if it is assumed to be zero then 
the expression for r direction velocity could be used to reduce the expression for the 
velocities in the 6 and 9 directions. The local pressure on the sphere surface is found from 
the Bernoulli Equation, and then integrated numerically around the entire surface to 
determine the net force on the sphere. 

Alternatively, the lift force could also be computed by using the equation given by 

Miloh, 

mm 

F^-fctPjrJl^X n(n+2) 2 A 11 A 11+1 «.2iip I rJu 2 f (56) 

n*l 

where f is the lifting force coefficient and is defined as 
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( 57 ) 


f= X n(n+2) 2 \ A n+ i 
n»l 

The lifting coefficients for two identical touching spheres with different number of 
terms in the velocity potential as computed by using Eq. (57) are listed in Table 3 and 
plotted in Fig. 14. 

Table 3. Lifting coefficients f for touching spheres with various 
number of terms n in the velocity potential 


N 

f 

N 

f 

10 

.4373060 

20 

.4921019 

30 

.5143508 

60 

.5404147 

100 

.5527179 

200 

.5633335 

250 

.5657013 

500 

.5706787 

1000 

.5738523 

2000 

.5755495 



Fig. 14 The lift coefficient (f) for two touching spheres or a sphere touching 
a plane against the numbers of terms in the velocity potential (N) 

Fig. 14 shows that as N increases f increases, but the increase levels off as n 

becomes large and f has a finite limit as N goes infinity. This limit can be estimated by 

using the data given in Table 3. It can be seen from Table 3 that each time n doubles the 

increase in f is less than 2/3 the increase it made previously, i.e. 

«4N) - «2N) < | («2N) -fTN)) 


( 58 ) 




where f(N) is the lifting coefficient f for a velocity potential having N terms. 

It is reasonable to assume that Eq. (58) will hold for N>2000, and the upper limit of f 


is then given as: 

f < f(2000)+ 1 (f(2000)-f(1000) ) + j (f(2000)-f(1000))+ (59) 

Upon solving the geometric series in Eq. (59), we have 

f<3K2000) - 2f(1000) (60) 

or 

f<5789421 (60') 

37 

which is approximately 77 . 

04 


It may be seen noted in Table 3 that the value of f for N=2000 is more than eleven 
per cent (11%) larger than that for N=30, while its value for N=2000 is less than .3% larger 
than that for N=1000. 

For two identical spheres with varying center distances, f computed by using Eq. 
(57) are plotted in Fig. 15, and some selected values and the corresponding center distance 
c are listed in Table 4 below . 

Table 4. Some selected lifting coefficients f for various center distance c 


c 

f 

c 

f 

2.000000 

£755125 

2.000001 

.5742419 

2.000002 

£732937 

2.000005 

.5713256 

2.000007 

J5703674 

2.000010 

.5691914 

2.000030 

1(642572 

2.000070 

.5585942 

2.000100 

.5555539 

2.000600 

.5347092 

2.001000 

.5206195 

2.005000 

.4685720 

2.010000 

.4345427 

2.060000 

.3187479 

2.100000 

.2526425 

2500000 

.0903779 

3.000000 

.0396850 

4.000000 

.0119794 

10.00000 

.0003003 

22.00000 

.0000128 
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Table 4 shows that when the center distance between two identical spheres is 10 
times their radius the interaction force between the two is almost negligible. Expressing 
this in another way, when the sphere center is 5 times its radius from a plane the sphere 
almost does not "see" the presence of the plane, and the interaction between them is 
negligible. 


VIL The Lift on a Hemisphere with the Base 
on an infinite plane parallel to the fluid velocity 


The velocity potential for this case is the same as a fixed sphere in an infinite ideal 
fluid flow field: 

U r b sinG 

= U r sinG cos<p + = — coso (61) 

2 r* 

The tangential velocities in the 0 and 9 direction can be obtained by 

differentiating the velocity potential as given by Eq. (61) with respect to 8 and 9: 

v 0 1 r-r = r lif ' r~r = UcosG CO89 + ^ COS0 CO89 = | Ucos6 cos<p (62) 

b b 

, 1 34> . ... . U . 3 .. . 

V« b = rsine 09 l r=T b =- Usin< P- 2'* in( P * (63) 

The fluid pressure on the surface of the hemisphere in the absence of earth gravity 
or when gravity is parallel to the plane to which the hemisphere base is attached is 
determined from the Bernoulli equation as: 

PI- ■- VV - p 0 -2 ' + 1 <») 

D Do 

* P 0 - 1 Pj U 2 (cos^ cos 2 9 ♦ sin 2 9 ) (64 ) 

Integrating the pressure distribution across the entire surface of the hemisphere 
and noting that the differential force is in the opposite direction to the surface normal, the 
Z-direction component of the force on the hemisphere is: 
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2 2x 


F t = - J" PI r _ T cos0 dA= it r 2 P Q + 1 Pj U 2 J J (cos 2 © cos 2 ^ +sin 2 9)r b cos0 sin0 dq> d0 


0 0 


*r, 


i 

' P o + I p i U2,tr b J (R 2 + DHdji.-ftr^+^xrJpjU 2 (65) 


Hie first term on the right in Eq. (65) represents the force when the pressure is 

uniform in the fluid and there is no flow, and the second term represents the lift force 

induced by the motion of the fluid around the hemisphere. The lift force coefficient in this 
27 

case is therefore gj. 

In the case where the direction of gravity is opposite to that of the base plane normal, 

the pressure on the surface of the hemisphere is: 

, 2 , 


Pl^= P© - 1 P| U 2 (cos 2 0 cos 2 <p + sin 2 9 ) - p ( gZ 


9 


P 0 - 1 Pj U 2 (cos 2 0 cos ^9 + sin z 9 ) - Pj g r b sinO 


.( 66 ) 


where P 0 is the pressure at the base plane, and the force component in the direction of the 


plane normal is: 

F z‘ * «&o + T r bPi& + § rtr b p ^ 2 


.(67) 


The first two terms on the right in (64) are the forces exerted on the hemisphere by 
the hydrostatic pressure and the third term on the right represents the dynamic force on the 
hemisphere as induced by the fluid flow. 

.When the gravity is in the same direction as the base plane normal, the pressure on 

the hemisphere surface is: 

Pl w * P 0 - 1 Pj U 2 (cob 2 0 cos 2 9 + «in 2 9 ) + p ( gZ 


-P 0 - g Pj U 2 (cos 2 © cos 2 9 + sin 2 9 ) + p, g r b sin© 

and the force on the hemisphere in the direction of the plane normal is: 

F i= ' xr b P o'T r b p i* + l ,lr bPi u2 


.( 68 ) 


.(69) 
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Similar to Eq. (67), the first two terms on the right in Eq. (69) are the hydrostatic force. 


and the third term is the dynamic force on the hemisphere. 

It can be seen from Eqs (65), (67) and (69) that the dynamic force or the lift on the 

hemisphere induced by the fluid motion is independent of the base plane orientation and the 

27 

lift coefficient for a hemisphere in a uniform flow field is 
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Fig. 2 RESIDUE VELOC 1 TY( IN PER CENT OF STREAM VELOCITY) AGAINST DEGREE ' 
FOR THE PARAMETER (NUMBER OF TERMS IN VELOCITY POTENTIAL) N-30 j 
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Fig. 2b RESIDUAL VE LOC I T Y ( i n percent of streom velocity) AGAINST DEGRE 
USING THE FIRST 30 TERMS OF A POTENTIAL HAVING 2000 TERMS 
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Fig. 7 RES I DUAL VELOC I T Y ( in percent of streom velocity) AGAINST DEGREE 
FOR THE PARAMETER (NUMBER OF TERMS IN VELOCITY POTENTIAL) N-2000 
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rig. 8 RESIDUAL VE LOC I T Y ( i n percent of streom velocity) AGAINST DEGRE 

FOR THE PARAMETER OF C/Rb-2 . 000001 ond N-2000 
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Fig. 10 RESIDUAL VEL0CITY(in percent of streom velocity) AGAINST DEGRE 

FOR THE PARAMETER OF C/Rb-2 . 000003 ond N-2200 
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Fig. 11 RESIDUAL VEL0CITY(in percent of streom velocity) AGAINST DEGREE 

FOR THE PARAMETER OF C/Rb-2 . 000005 and N-2500 


179 





Fig. 12 RESIDUAL VELOCITY(in percent of streom velocity) AGAINST DEGRE 

FOR THE PARAMETER OF C/Rb-2 . 000007 ond N-2300 
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Fig. 13 RESIDUAL VELOC I TY( i n percent of streom velocity) AGAINST DEGREE 

FOR THE PARAMETER OF C/Rb-2. 00001 and N-2200 
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A ppendix C. UNCERTAINTY ANALYSIS 


The uncertainties associated with the calculation of the gold film temperature, 
thermocouple temperature, system pressure, and total heat flux are estimated following the 
procedure of Kline et al. [114]. According to this procedure, the uncertainty u r in 
computing r is given by: 



(C.l) 


where r is a function of z\, Z2, .... z n and u Zj , u^, .... u^, are the uncertainties in each of 
these values, respectively. 

C . 1 Gold Film Temperature 

The CR7X was used for calibration of the gold film heater surfaces and for the 
imposed constant heat flux tests. The uncertainty in the heater surface resistance, which 
ultimately results in the heater surface temperature, was less for calibration of the heater 
surfaces than for actual power tests. From Equation (C.l) the uncertainty in the gold film 
resistance is expressed as: 


UR w 



(fuy 

v w 

2 


2 

f UR shl 

* 

V w 

VV w ) 

+ 

i v *j 

+ 

l R sh J 

) 


(C.2) 


Values of the quantities in equation from a representative calibration are: 


U V W = 0.6(1 V 

u V s h - °«^ v 

V w = 120mV 

V* - 50mV 

UR sh = 000 m 


R . = 3.760Q 

sh 



These values substituted into Equation (C.2) result in an uncertainty of 0.03%. For power 
runs, the above values are different and result in a greater uncertainty. Representative 
values for a power run are: 
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uy w = 0.05mV uvsh = 0.6jiV 

V w = 120mV = 50mV 

ur^ = 0.00000212 

R$h = 0.01631412 

with a resulting uncertainty in the heater surface resistance of 0. 1 1%. 

Calibration data with the largest scatter would be expected to provide a conservative 
estimate of the uncertainty in the slope. From the calibration data of surface Q-2, the 
corresponding temperature and resistance data with the largest deviation from the fit line are 
54.85°C at 2.8801 £2 and 83.59°C at 2.0171 £2. An estimate for the uncertainty in the 
slope of the linear relationship between the surface temperature and the resistance can be 
found by using Equation (C.l) together with the definition of this slope and the above 
values for temperatures and resistances, respectively. The uncertainty in the slope for these 
values is ± 0.23°C/£2. 

The uncertainty in the surface temperature measurement can now be estimated with a 
known estimate of the uncertainty in the heater surface resistance and of the slope in 
Equation (6). Carrying out the indicated operations of Equation (C.1) on the equation for 
the surface heater temperature, Equation (6), the following expression for uncertainty in the 
temperature measurement results: 

f 2 2 2M/2 

U T W - ((“T C ) + (“m) + (u(Rw-Rc)) J (C.3) 

Representative values for the quantities in equation are: 

U T C = 0.05°C ur w = 0.00312 

u m = 0.237c.£2 URc = 0.00312 
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These values result in an estimated uncertainty in the surface temperature measurement of 
about + 1 °C using the CR7X data acquisition system. 


C.2 Thermocouple Temperature 

The thermocouple spool from which all thermocouples used here were fabricated was 
calibrated against a platinum resistance thermometer which had an uncertainty in 
temperature of ± 0.00 1°C, and a calibration equation for the spool was determined from a 
fourth order polynomial least square fit of the measured thermocouple voltage and the 
determined platinum resistance thermometer temperature. In terms of the uncertainty in the 
CR7X voltage measurement, which was ± 0.6 pV, the uncertainty in the temperature 
determination is estimated to be better than ± 0.05°C (± 0.1°F). 

C.3 System Pressure 


There was little discrepancy, on the order of ± 0.005 psi between the Ruska pressure 
determination, with any uncertainty in pressure of ± 0.0003 psi and the calculated pressure 
resulting from the calibration equations. The uncertainty in pressure resulting from the use 
of these equations was then taken as ± 0.005 psi. From Equation C.1: 


up 

P 



(C.4) 


where P is the calculated pressure from the calibration equation and V is the pressure 
transducer voltage signal. For a typical run the following values were obtained: 
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